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Abstract 


To  date,  much  effort  has  been  spent  on  improving  the  detection  of  space  objects; 
however,  the  ability  to  track  and  catalogue  space  objects  remains  only  as  good  as  the  ability 
to  determine  the  object’s  position  and  angular  rate.  The  foundation  of  space  situational 
awareness  (SSA)  is  the  ability  to  detect  a  space  object  and  to  predict  its  location  in 
the  future.  In  order  to  accomplish  SSA  for  Geosynchronous  Earth  Orbit  (GEO)  space 
objects,  the  Defense  Advanced  Research  Projects  Agency  (DARPA)  developed  the  Space 
Surveillance  Telescope  (SST)  to  enable  ground-based,  broad-area  search,  detection  and 
tracking  of  small  GEO  objects  in  space.  In  general,  the  SST  tracks  along  the  sky  at  the 
sidereal  rate  where  the  stars  will  appear  to  be  stationary,  while  objects  within  space  not 
moving  at  this  rate  will  appear  to  move  across  the  focal  array  from  data  frame  to  data  frame. 
As  the  time  between  each  frame  is  known  and  the  position  of  the  detected  object  can  be 
determined  in  the  focal  array,  it  is  possible  to  measure  the  angular  position  and  angular 
rate  of  a  detected  object.  The  two  main  types  of  detection  algorithms  used  in  this  thesis 
are  the  binary  hypothesis  test  (BHT)  and  the  multi-hypothesis  test  (MHT),  which  both  rely 
on  a  log-likelihood  ratio  (LLR);  however,  the  MHT  algorithm  adds  an  additional  step  to 
correlate  nine  system  point  spread  functions,  or  hypotheses,  with  the  image  data.  These 
detection  algorithms  lead  to  varying  degrees  of  accuracy  and  precision  in  determining  the 
position  and  angular  rate  for  a  detected  space  object.  The  research  within  this  thesis  shows 
that  the  MHT  algorithm  is  more  accurate  and  precise  than  the  BHT  algorithm. 
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SURVEILLANCE  TELESCOPE 


I.  Introduction 


1.1  Space  Situational  Awareness 

Since  the  launch  of  Sputnik  in  1957,  the  need  for  Space  Situational  Awareness  (SSA) 
has  become  an  ever  more  important  task  in  conducting  space  operations  [17].  SSA  is 
predominantly  concerned  with  characterizing  and  tracking  human-created  objects  in  Earth 
orbit.  In  its  most  basic  sense,  SSA  is  the  ability  to  detect  a  space  object  and  predict  its 
location  in  the  future;  this  data  is  known  as  the  metric  data  of  an  Earth  orbiting  object  [17]. 
In  addition  to  the  metric  data,  SSA  seeks  to  characterize  objects  in  space  and  the  space 
environment. 

In  order  to  perform  SSA  using  electro-optics  systems,  three  steps  must  occur  in 
sequential  order:  1)  the  detection  of  an  Earth  orbiting  object,  2)  the  determination  of  the 
object’s  angular  position  and  rate,  and  3)  the  steps  to  determine  the  object’s  orbit  (orbital 
determination)  from  the  calculated  angular  position  and  rate  data.  The  detection  of  an  Earth 
orbiting  object  has  traditionally  been  done  with  three  main  types  of  sensors:  ground-based 
radar,  ground-based  electro-optics  (EO),  and  space-based  sensors  [17].  In  order  to  improve 
SSA,  the  detection  of  an  object  must  be  accompanied  by  at  least  some  degree  of  position 
information  on  that  object. 

The  detection  of  Near  Earth  Asteroids  (NEAs)  is  another  important  activity  that 
ground-based  electro-optics  are  used  for.  While  not  technically  part  of  SSA,  the  detection 
techniques  used  to  detect  satellites  and  other  man-made  debris  for  SSA  are  nearly  identical. 
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1.2  Electo-Optics 

The  data  used  in  this  thesis  was  primarily  gathered  from  the  Space  Surveillance 
Telescope  (SST),  which  is  an  EO  system  developed  by  Defence  Advanced  Research 
Projects  Agency  (DARPA).  The  main  advantage  of  EO  systems  over  radar  systems  for 
SSA  is  their  ability  to  more  easily  detect  objects  around  36,000  km  often  associated 
with  geosynchronous  orbits  [17].  While  radars  have  the  ability  to  determine  range  and 
range  rate  information  which  EO  systems  do  not,  it  requires  significant  power  SSA  radars 
to  search  for  objects  above  5,000  km  [17].  Passive  EO  systems  collect  light  reflected 
to  Earth  and  are  not  limited  by  the  same  power  restrictions.  With  an  increasingly 
crowded  geosynchronous  orbit,  which  contains  many  military  satellites  that  provide  critical 
capabilities  to  warfighters,  protecting  them  from  collision  with  space  debris,  meteors  and 
microsatellites  has  become  a  top  priority  [19]. 

In  general,  EO  systems  work  by  collecting  and  focusing  electromagnetic  radiation 
in  the  light  spectrum  reflected  by  an  orbital  object  into  an  image.  This  means  that  SSA 
EO  systems  only  work  when  the  target  is  illuminated  by  the  sun  and  the  telescope  is 
in  darkness  [17].  This  constraint  primarily  affects  the  detection  aspects  of  SSA.  Most 
detection  algorithms  work  by  finding  contrasts  in  light  intensity  from  the  background  light 
levels.  This  is  why  daylight  detection  remains  difficult  as  the  contrast  between  daylight 
light  levels  and  objects  trying  to  be  observed  is  very  small.  Most  often,  observations  made 
by  EO  systems  are  only  made  during  the  local  nighttime. 

1.3  The  Space  Surveillance  Telescope 

Most  Earth-based  telescopes  are  unable  to  see  objects  such  as  micro-satellites  and 
space  debris  that  threaten  military  satellites.  The  SST  was  developed  by  DARPA  beginning 
in  2002  to  overcome  these  limitations  by  enabling  ground-based,  broad-area  search, 
detection  and  tracking  of  small  objects  in  deep  space  for  purposes  such  as  space  mission 
assurance  and  asteroid  detection  [2].  The  US  Congress  mandated  that  National  Aeronautics 
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and  Space  Administration  (NASA),  with  the  help  of  the  Department  of  Defense  (DoD),  to 
catalog  by  the  year  2020  at  least  90  percent  of  comets  and  asteroids  larger  than  140  meters 
that  are  in  Earth  orbit  or  come  within  close  proximity  of  Earth  [4].  The  SST  was  designed 
to  scan  space  to  detect  and  correlate  unknown  space  objects  rather  than  dwell  on  stellar 
objects  over  relatively  long  periods  of  time  [20].  The  SST  uses  a  3.5  meter  aperture  which 
allows  for  significantly  better  performance  for  detecting  dim  objects  in  space  over  other 
current  SSA  optics.  Additionally,  the  large  aperture  allows  for  the  SST  to  have  a  wide  field 
of  view  allowing  the  SST  to  provide  full  sky  coverage  with  frequent  revisit  times  in  the 
hopes  of  detecting  and  cataloging  Earth  orbiting  objects  and  NEAs  within  its  field-of-view. 
In  general,  the  SST  tracks  along  the  sky  at  the  sidereal  rate;  this  means  stars  appear  to  be 
stationary,  while  objects  within  the  solar  system  appear  to  move  across  the  Charge-Coupled 
Device  (CCD)  array  from  frame  to  frame.  As  the  time  between  each  frame  is  known  and 
the  position  can  be  determined  in  the  CCD  array,  it  is  possible  to  measure  the  angular 
position  and  angular  velocity  of  a  detected  object.  Of  note,  the  data  used  in  this  thesis  was 
gathered  as  the  SST  was  in  a  test  mode  which  tracked  a  communications  satellite,  ANIK- 
Fl,  in  geostationary  orbit.  In  this  case,  the  stars  appeared  to  move  across  the  sky  at  the 
sidereal  rate  [20] . 

Many  of  the  current  SSA  telescopes  have  difficulty  finding  and  tracking  small  objects 
across  wide  tracks  of  sky  [7].  The  SST  was  designed  with  these  two  limitations  in  mind  to 
look  primarily  for  objects  in  geosynchronous  orbit  and  to  provide  a  much  wider  angle  (6 
square  degrees)  field-of-view  from  which  to  view  the  wide  tracks  of  the  sky.  The  SST  is  a 
3-mirror  Mersenne-Schmidt  type  telescope  with  a  3.5  m  primary  mirror  with  a  short  focal 
ratio  of  F/1.0  [19].  Additionally,  there  is  a  1.75  m  secondary  mirror;  this  secondary  mirror 
acts  like  a  stop  making  the  aperture  of  SST  appear  donut  shaped  with  an  outer  diameter  of 
3.5  m  and  in  inner  diameter  of  1.75  m.  Figure  1.1  shows  the  layout  of  SST’s  3  mirrors.  The 
compact  layout  of  the  SST  ensures  a  highly  agile  pointing  system. 
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Primary  Mirror 


Path  of  incoming  light 


'Tertiary  Mirror 


Secondary  Mirror 


Figure  1.1:  Optical  layout  of  the  SST  [19]. 


The  SST  also  has  a  6  square  degree  curved  CCD  focal  array  to  match  the  telescope’s 
inherently  curved  focal  surface  [19].  As  the  telescope  has  an  extremely  short  focal  ratio, 
F/1.0,  with  a  relatively  large  field-of-view,  the  curved  CCD  focal  array  allows  for  light 
entering  at  any  angle  within  the  6  square  degrees  FOV  to  be  focused  properly  at  the  3.5m 
focal  length  of  the  telescope.  The  6  square  degree  FOV  is  composed  of  3  x  2  degrees  FOV; 
to  capture  this  FOV,  the  curved  CCD  focal  array  is  12288  x  8192  pixels  [19].  Each  pixel  is 
a  square  \5fum;  thus  the  CCD  array  is  approximately  18.43cm  x  12.29cm. 

Note,  the  terms  pixel,  CCD,  and  photodetector  are  used  interchangeably  throughout 
this  thesis  and  refer  to  the  same  thing.  Additionally,  the  terms  intensity  and  irradiance 
are  also  used  interchangeably  throughout  this  thesis  and  are  considered  equivalent  for  the 
purposes  of  this  research. 

1.4  Detectors 

There  are  two  main  types  of  detection  algorithms  used  in  this  thesis.  The  first  type 
uses  a  Binary  Hypothesis  Test  (BHT),  which  computes  a  Log-Likelihood  Equation  (LLE) 
to  form  a  ratio,  often  referred  to  as  the  Log-Likelihood  Ratio  (LLR),  for  each  pixel  in  an 
image  [20].  The  LLR  is  the  ratio  of  an  object’s  likelihood  of  being  in  a  pixel,  denoted  as 
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hypothesis  Hu  to  the  likelihood  that  an  object  is  not  in  a  pixel,  denoted  as  hypothesis  H0. 
Each  hypothesis  is  given  equal  likelihood  of  occurring.  If  the  LLR  is  above  a  specified 
threshold  r,  then  H\  is  considered  to  be  more  likely  for  a  particular  pixel,  and  vice-versa,  if 
the  LLR  is  below  r,  then  H0  is  considered  to  be  more  likely.  The  probability  of  detection, 
Pd,  forms  an  inverse  relationship  with  the  probability  of  false  alarm,  P/a,  in  that  if  the 
threshold  value  that  determines  if  H{  or  H{)  is  more  likely  is  set  too  low,  then  the  Pja 
increases.  Lor  this  reason,  r  is  often  set  higher  than  1  in  order  to  reduce  the  Pfa  to  6-sigma. 
Currently,  the  BHT  is  the  detection  method  being  used  on  the  SST.  Section  3.3.1  of  this 
thesis  describes  the  BHT  method  in  greater  detail. 

The  second  type  of  detection  algorithm  uses  a  Multi-Hypothesis  Test  (MHT).  The 
MHT  is  identical  to  the  BHT  method  in  that  the  same  LLR  array  of  values  for  a  given  image 
are  computed;  however,  the  MHT  adds  an  additional  step  to  correlate  nine  hypotheses  with 
the  image  data  [20].  Each  hypothesis  is  a  sampled  Point  Spread  Lunction  (PSL).  The  use  of 
nine  PSPs  arise  from  the  fact  that  an  object  being  observed  is  not  always  focused  perfectly 
in  the  center  of  a  pixel.  If  an  observed  point  source  of  light  is  focused  to  the  middle  of 
two  pixels,  the  intensity  profile  will  be  integrated  over  both  pixel  areas  and  the  resultant 
PSP  will  be  appear  different  from  the  PSP  that  is  formed  if  the  observed  point  source  is 
focused  and  integrated  over  the  area  of  only  one  pixel.  The  nine  PSPs  are  generated  from 
the  various  combinations  of  shifting  the  observed  point  source  up,  down,  left,  and  right  by 
half  a  pixel  from  the  center. 

In  the  MHT,  the  nine  hypotheses  are  then  correlated  with  the  LLR  values.  The  nine 
generated  hypotheses  provide  a  more  accurate  representation  of  an  object  being  observed 
by  an  EO  system  and  can  thus  provide  a  more  confident  assertion  as  to  whether  an  object 
is  detected  or  not.  Additionally,  the  hypothesis  with  the  strongest  correlation  values  can  be 
assumed  to  represent  the  actual  observed  position  most  accurately.  This  additional  piece  of 
position  data  is  a  key  piece  of  information  that  the  MHT  provides  and  can  be  used  to  more 
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accurately  determine  the  angular  position  of  an  observed  object.  Section  3. 3. 2. 2  also  refers 
Multi-Hypothesis  Test  with  Outlier  Removal  (MHTOR),  which  contains  a  slightly  different 
method  from  MHT  in  computing  the  background  noise  levels  that  shows  some  additional 
improvement  over  the  more  traditional  MHT.  Section  3.3.2  contains  more  details  on  the 
MHT  methods.  Section  3.2  describes  the  optical  model  of  the  SST  and  the  generation  of 
the  nine  PSFs  in  more  detail. 

1.5  Determining  Position 

Once  an  object  is  detected,  then  step  2  of  the  SSA  process  can  occur,  which  is  to 
determine  the  angular  position  and  rate  of  a  detected  object.  The  algorithm  currently  used 
on  the  SST  to  determine  the  position  of  the  detected  object  is  to  centroid  the  binary  map 
of  ones  and  zeros  determined  by  the  BHT.  This  thesis  also  examines  two  other  options  for 
determining  the  position  of  the  detected  object.  The  first  alternate  option  is  to  use  a  Center 
of  Intensity  (Col)  algorithm.  This  option  uses  the  binary  map  of  the  image  determined  by 
the  BHT  and  maps  the  detection  points  back  to  the  intensity  values  of  the  original  image.  A 
centroiding  operation  is  then  performed  on  the  intensity  values  to  determine  the  Col.  This 
approach  has  shown  some  improvement  in  determining  position  over  centroiding  the  binary 
map  determined  by  the  BHT.  The  approach  taken  with  MHT  methods  is  to  centroid  the 
Signal  to  Noise  Ratio  (SNR)  values.  Only  four  of  the  hypotheses  from  the  MHT  were 
used  when  scanning  through  an  image  as  some  of  the  hypotheses  are  redundant.  The 
corresponding  correlation  values  for  the  four  hypotheses  were  then  input  into  a  matrix 
that  was  double  the  length  and  width  of  the  original  image.  The  centroiding  operation 
was  then  performed  on  this  large  matrix  which  has  all  the  correlation  values  from  the  four 
hypotheses.  The  results  from  this  last  method  seemed  to  provide  the  most  accurate  position 
information.  Section  3.4  discusses  the  position  algorithms  in  more  detail. 
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1.6  Measuring  Angular  Rate 

Once  the  position  data  has  been  determined,  it  is  possible  to  determine  the  angular 
rate  of  a  detected  object.  This  is  done  by  measuring  the  amount  of  change  in  an  object’s 
position  from  one  frame  to  the  next  which  is  measured  in  A  pixels  per  frame.  Using  the 
focal  length  of  the  optics  and  size  of  each  pixel,  it  is  possible  to  determine  what  each  pixel 
represents  in  radians.  This  ratio  is  then  used  to  convert  A  pixels  per  frame  into  radians  per 
frame.  As  the  time  between  each  frame  is  known,  the  angular  rate  in  radians  per  second 
can  be  determined.  Section  3.5  describes  this  process  of  determining  the  angular  rate  of  a 
detected  object  in  more  detail. 

This  thesis  primarily  uses  the  angular  rate  of  stars,  which  have  a  well  known  angular 
rate.  Stars  appear  to  orbit  the  Earth  at  the  sidereal  rate.  In  truth,  even  the  closest  stars 
are  over  at  least  4.4  light  years  away  and  their  position  does  not  change  with  respect 
to  Earth’s  rotation.  However,  as  the  Earth  rotates  once  per  sidereal  day,  the  stars  from 
a  topocentric  observation  point,  appear  to  orbit  the  Earth  once  per  sidereal  day  or  2 n 
radians  per  23.9344696  hours.  This  angular  rate  is  well  known  to  a  high  degree  of 
precision  and  accuracy  and  can  then  be  considered  as  the  truth  data  by  which  angular  rates 
calculated  as  part  of  this  research  can  be  compared  against.  The  least  squares  method  was 
used  to  determine  a  mean  and  variance  for  each  comparison.  Section  3.6  describes  the 
determination  of  the  mean  and  variance  using  the  least  squares  method  in  more  detail. 

1.7  Errors  in  Measuring  Position 

Additionally,  EO  systems  are  highly  subject  to  atmospheric  conditions  that  introduce 
phase  aberrations  to  the  propagated  wavefront.  Specifically,  the  tilt  phase  error  is  the 
strongest  phase  error  introduced  by  the  atmosphere  and  has  the  greatest  affect  on  the 
observed  object’s  pixel  position  [15].  The  determination  of  position  underlies  the  ability  to 
determine  angular  rate  and  must  be  considered.  Subsection  3.7.1  of  this  thesis  discusses  in 
more  detail  the  effects  of  the  tilt  phase  error  has  on  the  position  data  gathered  by  the  SST. 
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1.8  Problem  Statement 


In  a  purely  deterministic  system,  the  current  orbital  dynamics  of  an  observed  space 
object  are  presumed  to  perfectly  known  [18].  Based  on  these  orbital  dynamics,  it  is  possible 
to  know  exactly  where  an  object  will  be  in  the  future.  Under  this  ideology,  measurements 
of  where  a  space  object  is  are  presumed  to  be  perfect.  However,  atmospheric  turbulence 
and  imperfect  optics  cause  errors  in  the  observed  position  of  an  object.  With  atmospheric 
turbulence  being  a  random  phenomenon  and  the  inexact  measurement  of  optical  aberrations 
in  optics  system,  an  observation  uncertainty  when  measuring  the  position  of  an  observed 
object  must  be  taken  into  account.  Based  on  the  observation  uncertainty,  an  observed 
position  becomes  a  probability  that  the  object  is  at  a  position.  This  observation  uncertainty 
propagates  when  trying  to  predict  the  future  location  or  track  of  an  object  [18].  Therefore, 
when  measuring  the  position  of  an  object,  the  precision  and  accuracy  of  the  measurement 
are  very  important  parameters  for  determining  the  position  of  an  object  in  the  future.  In 
order  to  improve  SSA,  it  is  necessary  to  reduce  the  observation  uncertainty  by  improving 
the  accuracy  and  precision  of  the  measurement. 

1.9  Research  Objectives  and  Focus 

The  primary  research  objective  of  this  research  is  to  determine  which  detector  from 
both  the  point  detector  BHT  and  MHT  in  combination  with  a  centroiding  operation  are 
better  at  determining  the  position  of  a  detected  space  object.  While  research  has  been  done 
in  [20]  to  show  that  the  MHT  is  a  better  detector,  no  research  has  been  done  to  determine 
which  detector  is  better  able  to  provide  single  frame  position  information  from  which 
angular  rate  can  be  calculated  when  looking  at  changes  in  position  over  multiple  frames. 
The  research  will  focus  on  determining  the  angular  rate  of  a  star.  The  primary  reason  for 
determining  the  angular  rate  of  a  star  is  because  a  star  has  a  very  precise  and  known  angular 
rate  providing  accurate  truth  data.  As  each  algorithm  will  measure  an  angular  rate,  it  can 
be  compared  to  the  known  angular  rate  to  determine  which  algorithm  is  most  accurate. 
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As  each  algorithm  has  some  measured  error  from  the  actual  position,  the  determination  of 
which  algorithm  will  be  done  by  statistical  methods  -  specifically,  the  mean  and  standard 
deviation  will  be  determined  by  comparing  the  measured  data  against  the  truth  data.  Once 
a  particular  algorithm  is  determined  to  be  the  most  accurate  in  determining  the  angular  rate 
of  a  star,  then  it  can  then  be  used  to  determine  the  angular  rate  of  orbiting  objects  that  do 
not  have  a  known  angular  rate. 

1.10  Assumptions  and  Limitations 

The  main  limitation  for  the  algorithms  developed  in  this  research  are  that  they  were 
developed  by  windowing  the  data  from  the  SST  to  only  a  single  space  object  per  image 
frame.  As  the  centroid  operation  used  in  this  research  to  determine  position  would  be 
inaccurate  if  two  or  more  objects  were  detected  in  a  particular  window,  it  is  necessary  to 
window  the  detected  object  to  a  single  object.  There  are  clustering  operations  that  can 
separate  multiple  objects  and  thus  determine  position  and  angular  rates  of  multiple  objects; 
however,  it  is  outside  the  scope  of  this  research. 

Additionally,  as  the  algorithms  for  determining  position  were  only  verified  using  the 
angular  movement  of  stars,  which  have  a  known  angular  rate,  it  is  assumed  that  the  position 
algorithms  as  proposed  in  this  research  will  also  function  properly  in  determining  the 
angular  rate  of  an  object  that  has  an  unknown  angular  rate. 

1.11  Implications 

The  SST  uses  the  pixel  location  of  three  positive  detection  decisions  to  build  an  orbital 
track  of  the  object  which  can  then  be  used  to  predict  the  future  position  of  the  object  [20]. 
The  SST  currently  uses  the  integer  pixel  location  of  points  detected  using  the  BHT.  This 
approach  does  not  have  sub-pixel  location  from  which  to  create  the  track  of  the  detected 
object;  this  can  lead  to  large  errors  in  determining  the  track  of  the  object.  The  improvements 
in  the  detection  capability  using  the  MHT  have  been  documented  in  [20];  however,  this 
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detection  technique  also  allows  for  sub-pixel  location  data.  With  a  more  accurate  position 
algorithm,  as  proposed  in  this  research,  it  is  possible  to  reduce  the  position  error  of  each 
observation.  This  in  turn,  reduces  the  covariance  matrix  and  thus  better  estimates  of  the 
track  can  be  made. 

1.12  Preview 

This  research  aims  to  demonstrate  that  the  MHT  method  provides  better  position 
information  on  a  detected  space  object.  The  position  information  is  used  to  calculate 
the  angular  rate  of  the  detected  object;  thus,  the  better  position  information  will  lead  to 
determining  better  tracks  of  space  objects.  This  research  will  add  to  the  research  that  has 
already  been  done  in  [20]  to  show  that  the  MHT  is  a  superior  detector  over  the  currently 
used  BHT  on  the  SST  by  demonstrating  that  the  MHT  can  be  used  to  provide  better  position 
information  as  well.  Chapter  2  describes  the  process  of  how  the  SST  is  modeled,  the 
setup  of  both  the  BHT  and  MHT  methods  to  process  data  from  the  SST,  the  description 
and  setup  of  the  position  algorithms  used  in  this  research,  the  method  for  determining 
angular  rate  from  the  position  data,  and  the  method  for  determining  which  set  of  algorithms 
best  determined  the  angular  rate  from  a  statistical  perspective.  Chapter  3  will  interpret 
the  results  of  the  analysis  and  tests  run  using  the  position  algorithms  on  the  SST  data. 
Chapter  4  discusses  the  validity  of  the  proposed  position  and  angular  rate  algorithm  based 
on  comparison  to  the  current  method  being  employed  on  the  SST. 
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II.  Literature  Review 


2.1  Chapter  Overview 

This  section  of  the  paper  provides  a  survey  of  some  of  the  relevant  literature  that  has 
been  written  concerning  the  topics  of  1)  detecting  whether  an  object  is  present  or  not  in 
an  image  and  2)  determining  the  error  associated  with  making  dynamics  measurements 
using  topocentric  measurements  of  Earth  orbiting  objects.  Topic  1)  is  written  about 
by  Zingarelli,  Pearce,  Lambour,  Blake,  Peterson,  and  Cain  in  Improving  the  Space 
Sun’eillance  Telescope’s  Performance  using  Multi-hypothesis  Testing.  Topic  2)  is  written 
by  Maruskin,  Scheeres,  and  Alfriend  in  Correlation  of  Optical  Observations  of  Objects  in 
Earth  Orbit. 

2.2  Detection  Techniques 

The  Zingarelli,  et  al.  paper  describes  the  purpose  of  the  SST,  the  BHT  detection 
algorithm,  and  the  MHT  detection  algorithm  and  is  foundational  to  much  of  the  research 
done  in  this  paper.  The  point  of  the  experiment  run  in  the  Zingarelli,  et  al.  paper  is  to 
determine  which  type  of  algorithm  from  the  BHT  or  the  MHT  is  best  at  detecting  dim 
unresolvable  objects  in  space  on  a  single  frame  basis.  In  order  to  run  the  experiment  in 
the  Zingarelli,  et  al.  paper,  a  satellite,  ANIK-F1,  in  geosynchronous  orbit  that  is  gradually 
going  into  the  eclipse  behind  the  Earth  was  chosen.  This  allows  the  unresolvable  satellite 
to  experience  a  near-continual  decrease  of  solar  illumination,  thus  providing  a  continuum 
of  irradiance  values  over  which  to  test  the  performance  of  the  BHT  and  the  MHT.  The 
results  of  the  experiment  run  in  the  Zingarelli,  et  al.  paper  show  that  the  MHT  method 
improved  the  probability  of  detection  Pd  by  as  much  as  50  percent  over  the  point  detector 
BHT,  which  is  detection  algorithm  currently  run  on  the  SST. 
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2.2.1  Relevant  Research:  Detection  Techniques. 

The  SST  uses  at  a  minimum,  three  sequential  detections  to  determine  a  track  which 
can  be  used  to  predict  the  future  position  of  an  object.  The  track  of  an  object  is  based 
on  the  position  of  the  object  at  each  detection.  Therefore,  providing  the  most  accurate 
position  data  becomes  a  critical  piece  of  information  in  being  able  to  predict  the  future 
position  of  the  detected  space  object.  As  will  be  seen  in  this  research,  the  additional  sub¬ 
pixel  information  provided  by  the  MHT  will  provide  a  significant  advantage  in  determining 
the  pixel  position  information  over  the  BHT.  This  information  can  then  be  used  to  measure 
the  angular  position  and  angular  rate  of  an  object  which  is  the  primary  purpose  of  this 
paper.  Much  of  the  work  in  the  Zingarelli,  et  al.  paper  was  primarily  aimed  at  increasing 
the  detection  capability  of  the  SST  by  using  the  MHT  whereas  the  primary  purpose  of  this 
research  shows  that  while  the  MHT  increases  detection  capabilities,  it  also  provides  better 
position  information  of  the  detected  space  object  from  which  better  predictions  of  angular 
position  and  angular  rate  can  be  made. 

2.2.2  The  BHT  Detector. 

Currently,  the  SST  uses  a  BHT  algorithm  to  detect  space  objects  in  a  single  image 
[20].  This  approach  is  called  the  baseline  point  detector;  the  other  BHT  analyzed  in  the 
Zingarelli,  et  al.  paper  is  the  matched-filter  technique  (i.e.  correlation  detector).  Both  the 
point  detector  test  and  correlation  detector  are  derived  in  terms  of  SNR.  The  correlation 
detector  was  not  analyzed  in  this  thesis;  while  the  correlation  detector  is  shown  to  increase 
the  probability  of  detection,  it  does  not  provide  sub-pixel  location  information.  As  the 
primary  purpose  of  this  thesis  is  to  measure  the  angular  rate  of  a  detected  space  object,  the 
correlation  detector  is  not  analyzed. 

The  BHT  algorithm  works  by  setting  two  hypotheses  for  a  given  pixel  in  an  image. 
The  null  hypothesis  Hq  is  the  hypothesis  that  a  space  object  is  not  in  a  pixel.  The  alternate 
hypothesis  H\  is  that  a  space  object  is  in  a  pixel.  If  the  H\  hypothesis  is  more  likely  for  three 
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successive  frames,  then  an  image  is  considered  detected.  Additionally,  the  pixel  location 
of  the  three  positive  detection  decisions  is  used  to  build  an  orbital  track  of  the  object  which 
can  then  be  used  to  predict  the  future  position  of  the  object.  The  use  of  three  successive 
positive  detections  helps  eliminate  non-Gaussian  noise  sources  such  as  cosmic  rays  and 
data  communication  errors  arising  from  the  photo-detector,  thus  reducing  the  false  alarm 
rate  when  detecting  an  object. 

The  baseline  point  detector  used  in  the  SST  can  be  derived  from  a  Likelihood  Ratio 
Test  (LRT)  [20].  The  BHT  is  derived  in  terms  of  two  hypotheses.  The  hypothesis  H\  is  a 
likelihood  that  an  object  is  present  in  the  pixel  of  interest;  H0  is  the  hypothesis  that  an  object 
is  not  present  in  the  pixel  of  interest.  P(d(x,y)i(x,y)  e  [l,Md]\Hj)  is  the  joint  conditional 
probability  of  the  data  given  that  Hhi  e  {0, 1}  is  true.  Therefore,  in  the  BHT,  the  LRT 
becomes  a  ratio  of  the  hypotheses  probabilities.  In  order  to  minimize  the  probability  of 
false  alarm,  this  ratio  is  not  given  equal  weighting;  instead,  this  ratio  must  be  greater  than 
a  threshold  value  of  six  in  order  for  the  H\  hypothesis  to  be  considered  more  likely. 

2.2.3  Calculating  the  Probability  of  False  Alarm  for  the  BHT. 

The  threshold  value  inherently  sets  the  probability  of  false  alarm  PFA  [20].  While  a 
Poisson  distribution  is  a  more  accurate  representation  of  background  noise,  a  Gaussian 
distribution  can  be  used  to  approximate  the  image  noise.  In  the  case  of  the  SST,  the 
threshold  y  =  6  is  chosen.  Thus,  the  PFA  for  the  SST  is  defined  as  the  chance  that  a 
pixel  that  only  contains  background  light  (no  object)  will  produce  a  detector  output  that 
exceeds  the  threshold  value  of  6  [20].  When  an  object  is  not  present  in  a  pixel,  the  ratio 
produces  a  unit  variance,  zero  mean  Gaussian  random  variable;  therefore,  the  probability 
of  false  alarm  per  pixel  (ignoring  non-Gaussian  forms  of  noise)  is 

oo 

PFA=P(SNRbaseline>6\H0)=  f  —j=e^~dt  =  9.87c  -  010.  (2.1) 

J  yj2n 
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Table  2.1:  Alternative  Hypothesis  Sub-pixel  Shifts 


Alternative(a) 

Horizontal  Shi  ft  (o' A 

Vertical  Shifty,) 

1 

0 

0 

2 

0 

-15  /mi 

3 

0 

+  15  nm 

4 

+  15  urn 

0 

5 

-15  nm 

0 

6 

+  15  nm 

-15  nm 

7 

+  15  nm 

+  15  /urn 

8 

-15  nm 

-15  /im 

9 

-15  nm 

+  15  nm 

2.2.4  The  MHT  Detector. 

In  the  Zingarelli,  et  al.  paper,  another  matched  filter  method,  the  MHT,  is  proposed  for 
single  frame  detection.  The  primary  rationale  for  the  introduction  of  the  MHT  is  because 
the  image  of  a  space  object  does  not  always  get  focused  in  the  center  of  a  pixel.  As  the 
shape  of  a  sampled  PSF  changes  depending  on  where  the  object  is  imaged  on  the  CCD 
array,  the  correlation  with  a  system  PSF  with  the  image  data  may  not  correlate  strongly. 
In  order  to  account  for  the  possibility  that  the  space  object  may  be  focused  off-center  of 
a  pixel,  a  multiple-hypothesis  strategy  is  considered.  Each  hypothesis  is  a  sampled  PSF 
which  is  generated  by  shifting  the  PSF  by  plus  and  minus  a  half-pixel  (which  is  15 /urn  in 
the  SST  CCD  array)  horizontally  and  vertically  as  described  in  Table  (2.1);  the  system  PSF 
is  then  sampled  to  create  the  nine  hypotheses.  Each  hypothesis  //,  is  then  correlated  with 

the  data  in  each  pixel.  The  MHT  method  selects  from  nine  hypotheses  (H j  -  H9)  and  the 
null  hypothesis  H0  to  determine  whether  the  image  is  in  the  center,  a  corner,  or  a  side  of 
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a  pixel.  The  Multi-hypothesis  decision  to  select  one  hypothesis,  Hk,  over  another  one  of 
the  other  nine  hypotheses,  Hh  is  based  on  uniform  cost  and  equal  priors  using  conditional 
probabilities  [20]. 

2.2.5  Calculating  the  Probability  of  False  Alarm  for  the  MHT. 

While  the  goal  of  this  thesis  is  not  specifically  in  increasing  the  detection  of  space 
objects,  it  is  important  to  consider  the  probability  of  false  alarm  of  the  MHT.  Using  the 
Eq.  (2.1)  and  considering  each  hypothesis  to  be  mutually  exclusive  from  one  another  and 
also  assuming  that  the  result  of  each  individual  sub-MHT  to  be  statistically  independent  of 
each  other,  the  PFA  is  ~  4x9. SI e  -  10  =  3.94<?  -  9  [20].  The  estimated  PFA  is  higher  for  the 
MHT  than  for  the  BHT;  however,  the  PFA  can  be  reduced  by  raising  the  detection  threshold 
Jmht  =  6.2212,  which  yields  a  PFA  =  9.87<?  -  10. 

2.2.6  MHT  With  Outlier  Removal. 

The  Zingarelli,  et  al.  paper  also  proposes  a  separate  method  for  computing  the 
standard  deviation,  cr,  of  the  background  noise  as  listed  in  Eq.  (3.20).  This  method  rejects 
any  noise  sample  in  the  window  surrounding  the  pixel  to  be  tested  whose  values  are  outside 
those  predicted  by  Gaussian  statistics.  This  has  the  effect  of  eliminating  extremely  high 
noise  or  nearby  stars  in  the  calculation  of  background  noise. 

2.2. 7  Modeling  the  Point  Spread  Function  for  the  SST. 

In  the  Zingarelli,  et  al.  paper,  a  modeled  system  PSF  is  used  to  characterize  the  SST’s 
impulse  response,  hoptics(x,  y).  The  SST’s  optics  are  modeled  as  linear  and  shift  invariant. 
The  light  propagating  from  a  distant  point  source  such  as  a  star  or  Earth  orbiting  object  is 
assumed  to  be  temporally  incoherent  light.  The  Zingarelli,  et  al.  paper  describes  a  method 
for  generating  the  system  PSF  by  taking  the  inverse  Fourier  transform  of  the  product  of  the 
optics,  sampling,  and  atmospheric  transfer  functions.  This  method  is  also  utilized  in  this 
research. 
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2.2.8  Summary  and  Relation  to  this  Research. 

The  Zingarelli,  et  al.  paper  provides  the  method  for  implementing  the  BHT,  the  MHT, 
and  creating  a  system  PSF.  The  paper  also  states  that  the  MHT  improved  the  detection 
capability  of  the  SST  over  that  of  the  currently  employed  point  detector  test  while  also 
providing  sub-pixel  position  information.  This  research  utilizes  much  of  this  paper  to 
implement  the  BHT  and  MHT.  Additionally,  this  research  builds  upon  the  research  of 
this  paper  by  utilizing  the  sub-pixel  position  information  to  determine  the  angular  rate  of 
a  space  object.  By  utilizing  the  better  position  information  of  the  MHT,  this  research  also 
shows  that  the  track  information  that  the  SST  calculates  could  be  improved  by  utilizing  the 
MHT  method. 

2.3  Determining  the  Error  Associated  with  Making  Orbital  Dynamics  Measure¬ 
ments 

The  Maruskin,  et  al.  article  describes  a  method  for  determining  the  orbits  of  medium 
to  high  Earth  orbiting  objects  by  determining  the  intersection  of  two  orbital  tracks  from  two 
observation  points.  Due  to  the  large  uncertainty  associated  with  measuring  and  estimating 
range  and  range  rate  values  of  a  detected  object,  orbital  determinations  based  on  these 
estimates  will  likewise  also  have  large  uncertainties.  The  Maruskin,  et  al.  article  devises 
an  algorithm  which  does  not  rely  on  estimates  of  range  and  range  rate  to  perform  the  initial 
orbit  determination. 

From  a  mathematical  perspective,  once  the  angular  rates  and  accelerations  are 
measured  for  both  the  right  ascension,  a,  and  declination,  8,  angles,  it  is  possible  to 
determine  the  range,  p,  and  range  rate,  p  [11]. 

2.3.1  Determining  the  Error  in  Measuring  Angular  Position,  Rate,  and  Accelera¬ 
tion  and  the  Observation  Uncertainty. 

The  main  value  gained  from  the  Maruskin,  et  al.  article  is  in  determining  the  errors  in 
measuring  the  angular  position,  rate,  and  acceleration  of  a  space  object  due  to  observation 
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uncertainties.  Assuming  that  the  angular  motion  can  be  modeled  kinematically  and  that 
the  orbit  of  the  detected  space  object  is  deterministic  in  nature,  then  it  is  possible  to  fully 
determine  the  orbit  of  an  object  [11].  However,  from  a  practical  perspective,  there  is  a 
large  uncertainty  in  measuring  the  angular  position,  rate,  and  acceleration  of  an  orbital 
object  which  translates  into  a  large  uncertainty  in  estimating  the  range  and  range  rate  [11]. 
By  using  the  equations  described  in  the  Maruskin,  et  al.  article  and  further  elaborated 
on  in  section  (3.7.2),  it  is  possible  to  determine  the  error  associated  with  measuring 
angular  position,  rate,  and  acceleration  using  the  SST.  The  calculation  of  the  observation 
uncertainty  parameter  used  in  the  calculation  of  the  errors  in  the  angular  parameters  is 
elaborated  further  in  section  (3.7). 

2.3.2  Maruskin:  Summary /How  my  research  is  different. 

While  the  determination  of  the  orbit  of  a  detected  object  for  the  purposes  of  cataloging 
the  object  is  one  of  the  primary  goals  of  SSA,  the  method  of  correlating  multiple  optical 
observations  as  described  in  the  the  Maruskin,  et  al.  article  is  still  a  new  and  relatively 
untried  method.  Additionally,  the  method  described  in  the  article  relies  on  optical 
observations  that  are  several  minutes  apart  and  from  multiple  observation  points;  this 
is  a  departure  from  how  the  SST  operates  by  taking  multiple  sequential  measurements 
separated  by  as  little  as  0.675  seconds  [19].  For  these  reasons,  the  orbital  determination 
method  described  in  this  article  was  not  utilized  in  this  research.  However,  the  method 
for  calculating  the  error  in  the  angular  parameters  is  very  applicable  to  the  research  of  this 
paper. 
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III.  Methodology 


3.1  Chapter  Overview 

The  key  purpose  of  the  research  in  this  thesis  is  determining  the  angular  rate  of  a  space 
object,  thus  determining  the  most  accurate  pixel  position  of  an  object  detected  by  one  of  the 
matched- filter  tests  is  imperative.  The  SST  currently  uses  the  point  detector  BHT  to  detect 
a  space  object  [20].  As  the  BHT  name  indicates,  the  pixels  that  an  object  is  detected  in  are 
indicated  with  a  1,  while  the  pixels  that  no  object  is  detected  in  are  indicated  with  a  0,  thus 
creating  a  binary  map  of  the  sky.  If  a  single  space  object  detected  covers  more  than  one 
pixel,  then  a  centroid  operation  is  used  to  determine  the  position  of  that  object;  otherwise, 
the  integer  value  of  the  detected  pixel  is  the  position  of  the  object.  This  method  fails  to  take 
into  account  the  intensity  or  SNR  values  of  a  detected  object  therefore  introducing  errors 
into  the  calculation  of  position  of  a  detected  object.  Furthermore,  the  error  in  position 
measurement  degrades  the  accuracy  from  which  the  track  of  the  object  is  determined  from 
the  position  information  of  three  subsequent  positive  detections. 

Two  other  methods  were  analyzed  in  this  research  and  tested  with  data  from  the  SST 
to  determine  the  position  of  an  object;  these  included  performing  a  Col  algorithm  on  the 
BHT  and  using  the  SNR  values  from  MHT  and  centroiding  these  values.  The  Col  method 
maps  the  detection  points  from  the  BHT  method  back  to  the  original  data  which  contains 
the  intensity  values.  The  intensity  values  of  the  pixels  immediately  adjacent  to  the  detection 
point  were  then  centroided  to  determine  the  Col  and  position  of  the  detected  object.  The 
other  method  utilized  in  this  research  to  determine  the  position  of  a  detected  object  is 
based  on  the  MHT  method.  The  results  in  the  Zingarelli,  et  al.  article  show  that  the  MHT 
method  can  increase  the  probability  of  detection  of  a  space  object  over  the  BHT  while 
also  providing  sub-pixel  location  data.  In  the  MHT  method  for  determining  position,  the 
hypothesis  with  the  highest  SNR  is  selected  as  the  most  probable  hypothesis.  The  method 
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of  centering  the  SNR  values,  referred  from  this  point  onward  as  Center  of  SNR  (CoS), 
relies  on  taking  the  array  of  SNR  values  as  calculated  from  each  of  the  hypotheses  in  the 
MHT,  mapping  these  to  a  master  array  of  SNR  values  from  four  of  the  hypotheses,  and 
then  performing  a  centroid  operation  on  the  master  array  to  find  the  position. 

The  following  sections  in  this  chapter  explain:  1)  the  modeling  of  the  system  PSF 
for  use  in  the  MHT  method,  2)  the  setup  of  the  BHT  and  MHT  detectors,  3)  the  methods 
for  determining  position  of  a  detected  object  using  the  centroid  operation  on  the  binary, 
intensity,  and  SNR  values  as  determined  from  the  BHT  and  MHT  methods,  4)  the  method 
for  measuring  angular  rate,  5)  the  statistical  method  to  compare  each  of  the  angular  rate 
determination  methods  against  each  other,  and  6)  the  method  for  measuring  the  observation 
error  of  the  SST  and  the  errors  associated  in  measuring  the  angular  position  and  rate  of  a 
detected  object. 

3.2  Modeling  the  PSF  of  the  Space  Surveillance  Telescope 

To  create  the  simulated  PSFs,  the  SST  system  is  modeled  in  Matlab  by  using  the 
aberration  free  (diffraction-limited)  incoherent  Optical  Transfer  Function  (OTF)  of  the  SST. 
The  phase  aberrations,  which  are  departures  from  an  ideal  spherical  wavefront  form  due  to 
imperfections  in  the  optics  of  the  SST,  were  then  added  to  the  idealized  pupil  function  to 
create  the  generalized  pupil  function,  £P(x,y).  The  phase  aberrations  were  modeled  using 
Zernike  polynomials  and  were  provided  by  the  MIT  Lincoln  Labs.  The  resulting  intensity 
PSF  is  then  integrated  over  the  area  of  a  pixel  and  downsampled  to  the  SST’s  actual  pixel 
size.  This  last  step  creates  a  pixelated  version  of  the  PSF  known  as  the  sampled  PSF  in  this 
research. 

The  ability  to  accurately  model  the  PSF  of  the  SST  is  crucial  to  being  able  to  run  the 
MHT,  which  is  a  matched- filter  test  relying  on  the  cross-correlation  of  multiple  hypotheses. 
These  hypotheses  are  alternately  sampled  PSFs,  with  the  image  data  from  the  SST.  The 
underlying  assumption  for  the  MHT  to  be  considered  valid  is  that  most  objects  detected  by 
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the  SST  are  far  enough  away  and  small  enough  that  they  can  be  treated  like  a  point  source 
when  observed  from  the  SST.  If  this  assumption  is  valid,  then  the  data  gathered  from  the 
SST  can  be  treated  as  if  it  is  from  a  point  source;  therefore,  the  cross-correlation  of  the  data 
and  the  hypotheses  is  valid.  In  order  to  show  that  detected  objects  can  be  treated  as  a  point 
source,  a  simple  comparison  of  two  similar  triangles  as  in  Figure  3.1  is  used.  In  the  figure, 
3.5m  is  the  focal  length  of  the  SST,  30 /wi  is  the  effective  pixel  size,  and  36,000km  is  the 
approximate  distance  to  geosynchronous  orbit.  Therefore,  the  approximate  maximum  size 
that  a  space  object  at  geosynchronous  orbit  can  be  and  still  fit  within  one  pixel  is  308.6m. 
As  most  space  objects  that  SST  is  trying  to  detect  are  at  geosynchronous  orbit  are  smaller 
than  308.6m,  it  can  be  assumed  that  most  objects  can  be  treated  simply  as  a  point  source. 
While  some  asteroids  may  be  much  larger  than  308m,  they  are  at  a  much  higher  altitude 
than  geosynchronous  orbit,  and  thus  still  remain  within  the  requirements  for  being  treated 
as  a  point  source.  Without  this  assumption,  the  cross-correlation  of  an  imaged  space  object 


Figure  3.1:  Using  similar  triangles  to  calculate  maximum  size  of  a  space  object.  Note,  the 
figure  is  not  to  scale  and  is  drawn  for  clarification. 


and  the  SST  PSF  would  not  be  valid. 


20 


Now  that  it  can  be  assumed  that  a  detected  object  can  be  considered  to  be  a  point 
source,  the  determination  of  the  system  PSF  for  the  SST  can  be  considered.  In  the 
Zingarelli,  et  al.  paper,  a  modeled  PSF  is  used  for  the  characterizing  SST’s  impulse 
response,  hoptics(x,  y)  and  is  considered  the  optical  PSF.  The  SST’s  optics  are  modeled 
as  linear  and  shift  invariant.  The  light  propagating  from  a  distant  point  source  such  as  a  star 
or  Earth  orbiting  object  is  assumed  to  be  temporally  incoherent  light. 

Based  on  a  linear  and  shift- invariant  model  for  the  optical  system,  the  image 
irradiance,  i(x,  y),  in  the  CCD  array  plane  with  coordinates  (x,  y),  of  the  continuous  image  is 
the  convolution  of  a  point  source,  mathematically  expressed  as  a  delta  function,  6(x,y),  with 
a  telescope’s  impulse  response,  hoptics(x, y)  [20].  To  most  accurately  model  a  system’s  PSF 
it  is  necessary  to  take  into  account  the  wavefront  errors  introduced  by  the  optics  aberrations 
into  pupil  function,  (u,  v),  where  u  and  v  represent  coordinates  in  the  pupil  plane.  The 
wavefront  errors  are  represented  mathematically  as  the  first  N  Zernike  polynomials, 
where  N  is  the  maximum  number  of  polynomials  needed  for  a  given  telescope  [20]. 
The  wavefront  error  is  then  captured  by  scaling  the  Zernike  polynomials  with  a  Zemike 
coefficient,  Z2  -  ZN.  The  total  wavefront  error,  <5,  is  then 


S(u,v)  =  Z20  2(w,  v)  +  ...  +  ZN(pN(u,v).  (3.1) 

The  generalized  pupil  function  taking  into  account  telescope  aberrations  is  then  represented 
by 

P(u,v)  =  Jl(u,v)e[j£M\  (3.2) 


To  determine  the  incoherent  PSF,  the  field  at  the  pupil  is  then  propagated  using  a  Fraunhofer 
propagator  to  the  CCD  array,  the  square  of  the  field  is  then  taken  to  determine  the  irradiance 
[20].  The  irradiance  is  the  physically  measurable  attribute  of  an  optical  wave  field  [9].  The 
irradiance  pattern  is  the  PSF  of  the  system,  which  is  mathematically  represented  by 


hoptics(.Xi  y) 


l  l  j2n(xu+yv) 

I  I  P(u,v)e  v  dudv 


(3.3) 
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where  z  is  the  propagation  distance  and  A  is  the  mean  wavelength  of  the  light  wave  [9]. 
The  optical  transfer  function,  <Hoptics(fx,  fy),  is  then  computed  by  performing  a  Fourier 
transform,  T , 

U>Pacs(U  fy )  =  T  {hoptkJx,  >’)}  •  (3 .4) 

The  resulting  PSF  must  then  be  sampled  to  account  for  finite  square  pixels  of  the  SST, 
which  are  a  =  30 pm  [20].  Each  pixel  is  represented  by  a  rectangle  function.  This  allows 
the  irradiance  pattern,  which  describes  the  number  of  photons  arriving  at  a  location  (u,  v) 
in  the  CCD  array,  to  be  integrated  over  the  square  pixel.  The  total  power  s]3  incident  on 
the  surface  area  of  a  photodetector  is  a  function  of  the  responsivity  of  the  detector  and  the 
direction  of  power  flow  k  [9] .  This  is  represented  mathematically  as 

s]3  =  f  f p{u,v)—~dudv,  (3.5) 

J  J  \k\ 

A 

where  h  is  a  unit  vector  pointing  into  the  surface  of  the  detector  and  A  is  the  area  of 
the  photodetector.  When  k  is  nearly  normal  to  the  surface,  as  is  the  case  with  the  SST, 
T  is  simply  the  integral  of  the  power  density  p  over  the  detector  area  [9].  The  s]3  then 
represents  the  total  amount  of  photons  that  arrive  over  a  surface  area  A  over  a  given  amount 
of  time.  The  operation  can  be  represented  by  a  transfer  function,  fHPixei ,  which  is  the  Fourier 
transform,  T' ,  of  the  rectangle  function 

pixeiifx ,  fy)  =  T  {rect(ax,  ay)} .  (3 .6) 

Additionally,  atmospheric  turbulence  affects  the  irradiance  pattern.  As  the  integration 
time  used  on  the  SST  is  greater  than  25ms,  a  long-exposure  atmospheric  transfer  function 
is  an  acceptable  model  for  atmospheric  effects  [20] .  The  atmospheric  transfer  function  is 

KauniUfy)  =  exp  | -3 .44^  |  (3.7) 

where  the  radial  frequency  v  -  (  /)2  +  /'2)l/2  and  r0  is  the  atmospheric  coherence  diameter 
or  seeing  parameter  [8]. 
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The  overall  system  impulse  response,  hmodei(x,  y )  centered  on  a  pixel  is  then  the  inverse 


Fourier  transform  of  the  product  of  all  the  individual  transfer  functions: 

hmodel^X*  y)  —  1  [HopticsiU  fyYHptxel{U  fyYKumiU  fy)}  •  (3.8) 


A  Nyquist  sampled  PSF  would  not  change  the  shape  of  the  image  irradiance  pattern  as 
a  function  of  a  sub-pixel  shift,  A.v  [20].  However,  as  the  SST  is  undersampled  and  the  point 
source  is  not  always  focused  in  the  center  of  a  pixel,  the  irradiance  pattern  does  change  as  a 
function  of  where  the  point  source  is  focused.  To  model  these  different  irradiance  patterns, 
the  modeled  PSF  is  down-sampled  using  the  ratio  L  between  the  30 /urn  pixels  of  the  SST 
and  the  Nyquist  pixels  size  from  Eq.  (3.3).  The  shifted  and  down-sampled  PSF  is 


hsampim,  Ax,  n,  Ay) 


hmodei(x,  y)6(Lm  -  x  -  LAx,  Ln-y  -  LAyjdxdy, 


(3.9) 


and  the  sampled  irradiance  is 


isamp(m.  Ax ,  n.  Ay)  =  6hsamp(m ,  Ax,  n,  Ay)  +  B,  (3.10) 

where  B  is  the  background  light,  and  6  is  the  total  number  of  photons  emitted  from  the 
object  per  integration  time,  and  m  is  the  integer  value  for  each  pixel  location  in  the  CCD 
array  [20].  Using  the  shifts  as  described  in  Table  (2.1),  the  nine  different  hypotheses  used 
in  the  MHT  in  Eq.  (3.30)  were  determined. 

In  general,  the  method  for  generating  the  system  PSF  in  this  research  is  left  unchanged 
from  what  is  done  in  [20].  This  method  for  generating  the  system  PSF  is  based  on 
multiplying  the  optics,  atmospheric,  and  sampling  transfer  functions  together  to  form  the 
total  system  PSF;  this  is  summarized  in  the  following  equation: 

hmodel(x,y)  r-1  [H„PnCS(U  fyfHplxeifx,  fyfHaM,  fy)]  •  (3.11) 

Eqs.  (3.4),  (3.6),  and  (3.7)  are  used  to  create  the  total  system  impulse  response  hmodei(x,y); 
this  is  the  total  system  or  modeled  PSF. 
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The  pupil  function,  !P(w,  v),  used  in  Eq.  (3.3),  is  created  in  MATLAB  by  creating  a 
2400  x  2400  array  of  ones  and  zeros  to  represent  the  aperture  of  the  SST  where  the  ones 
represent  the  aperture  opening.  The  aperture  of  SST  is  donut  shaped  with  a  3.5m  outer 
diameter  and  a  1.75m  inner  stop;  Figure  (3.2)  depicts  the  aperture  of  the  SST  as  represented 
in  MATLAB. 


Figure  3.2:  The  aperture  of  SST  in  a  2400  x  2400  array. 


The  size  of  the  aperture  has  a  diameter  of  1200  samples,  which  is  half  of  the  area  the 
total  array.  The  sampling  size  is  set  at  2400  x  2400  samples  to  meet  Nyquist  sampling 
requirements  to  avoid  aliasing.  However,  the  Fourier  transform  of  a  circular  aperture,  as  is 
the  case  with  the  SST,  will  result  in  a  Bessel  function  J i  of  the  first  kind  represented  by 

T{circ{r)}  =  (3.12) 

P 

where  r  represents  the  radius  in  the  (x,  y)  plane  and  p  is  the  radius  in  the  two-dimensional 
spatial  frequency  plane  (w,  v)  [8].  That  being  said,  the  cutoff  frequency,  fc ,  is  set  by  the 
diameter,  D ,  of  the  aperture,  the  focal  length,  Zf,  and  the  mean  wavelength  A  of  the  light 
being  reflected  or  emitted  in  the  following  relationship  [9] : 

D 

fc  =  (3.13) 

Azf 
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According  to  Nyquist  sampling  theorem,  the  sampling  frequency  must  be  at  least  twice  that 
of  the  cutoff  frequency  (fs  >  2 fc).  Therefore,  the  sampling  frequency  fs  is: 


(3.14) 


To  satisfy  the  Nyquist  sampling  theorem,  a  2400  x  2400  sampling  array  is  selected,  which 
is  twice  the  aperture  diameter  of  1200  samples.  It  is  also  worth  noting  that  by  using  a 
A  =  500nm  and  given  that  Zf  =  D  =  3.5 m  for  the  SST,  the  sampling  period  Ts  is  0.25 fim. 
For  the  SST,  which  has  an  effective  pixel  size  of  30/rm,  this  means  that  there  are  120 
samples  per  pixel  per  a  given  coordinate  axis.  In  general  for  the  SST,  the  given  window 
that  the  background  noise,  B ,  is  calculated  from  is  about  a  20  x  20  pixel  window.  When 
considering  the  2400  x  2400  sampling  array  used  for  the  P(u,  v)  is  divided  by  120  in  either 
coordinate  direction,  it  translates  nicely  into  a  20  x  20  pixel  window.  This  was  purposely 
set  to  meet  the  requirements  of  Whittaker-Shannon  sampling  theorem,  which  implies  that 
the  exact  recovery  of  a  bandlimited  function  can  be  achieved  from  an  appropriately  spaced 
rectangular  array  of  its  sampled  values  [9].  The  two-dimensional  Fast  Fourier  transform 
was  then  performed  in  MATLAB  on  the  2400  x  2400  array  representing  the  pupil  function, 
rP(u,  v),  of  the  SST  to  determine  the  diffraction  limited  OTF,  Hoptics. 

Next,  in  order  to  create  an  accurate  model  of  the  system  PSF,  the  aberrations  or 
errors  in  the  wavefront  due  to  imperfections  in  the  optics  are  included.  An  aberration 
is  any  deviation  from  the  ideal,  diffraction-limited  imaging  performance  of  a  lens  [1]. 
Additionally,  the  Fraunhofer  approximations  most  significantly  impact  the  phase  factor 
of  the  spherical  wave  being  propagated  by  approximating  them  to  parabolic  wavefronts 
[9].  Mathematically,  lens  aberrations  can  be  described  by  Zemike  polynomials,  which 
help  express  the  exit-pupil  wavefront  of  the  light  being  propagated  [3].  According  to  the 
Woods,  et  al.  paper,  the  SST’s  Zemike  coefficients  for  the  first  twelve  Zemike  polynomials 
were  estimated  using  the  DONUT  algorithm,  which  generates  an  estimate  of  the  wavefront 
aberrations  by  directly  fitting  a  model  of  the  aberrated  image  to  an  observed  defocused  star 
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Table  3.1:  Mean  of  36  DONUT  Algorithms  Run  on  the  SST  [19]. 


Noll  index  (j) 

order  (n) 

frequency  (m) 

Classical  name 

Mean  Zernike  Coefficient 

1 

0 

0 

piston 

2.0877 

2 

1 

-1 

x-tilt 

-5.9517 

3 

1 

1 

y-tilt 

-5.2999 

4 

2 

-2 

defocus 

6.8849 

5 

2 

0 

y-astigmatism 

1.2598 

6 

2 

2 

x-astigmatism 

-0.2782 

7 

3 

-3 

y-coma 

0.2755 

8 

3 

-1 

x-coma 

-0.7279 

9 

3 

1 

y-trefoil 

0.3548 

10 

3 

3 

x-trefoil 

-0.4818 

11 

4 

-4 

spherical 

-0.1567 

[19].  The  DONUT  method  was  performed  on  thirty-six  defocused  stars;  the  mean  of  these 
36  tests  is  used  as  Zernike  coefficient.  Table  (3.1)  contains  the  mean  Zernike  coefficient 
for  each  Zernike  polynomial.  The  equation  to  model  the  Zernike  polynomials  is  based 


on  a  pyramid  where  the  rows,  k,  are  the  order  and  the  columns,  rj,  are  the  frequency.  The 
polynomials  are  defined  as 


where 


Zeven.  =  Vk  +  1  Kir)  \f2cos(i]0)  i]  f  0 
Zocidj  =  Va-+  1  R'l(r)  yflsinind)  //  f  0 
Zj  =  ^ Ik  +  \R°K(r),  rj  =  0 


(k-tD/2 

*?«=  z 


(-!)'(* -J)! 


,.k-2s 


s= 0 


5!  [(a-  +  if)/ 2  -  5]!  [(/c  -  rj)/ 2  -  .?]! 


(3.15) 


(3.16) 
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Figure  (3.3)  represents  the  phase  error  listed  in  Eq.  (3.1)  pictorially.  Figure  (3.4)  shows  the 
non-downsampled  PSF  from  the  SST;  this  is  the  irradiance  pattern  of  Eq.  (3.3),  which  only 
takes  into  account  the  pupil  function  and  the  phase  errors  of  the  SST. 


x-pixel 

Figure  3.3:  The  optics  phase  error  (perror  as  modeled  by  Zemike  polynomials. 


Next,  the  atmospheric  transfer  function,  Hatm(fx,fy),  is  modeled  for  the  data  from 
the  SST  using  Eq.  (3.7).  The  seeing  parameter,  r0,  for  the  data  used  in  this  research  is 
determined  using  a  phase  retrieval  method  described  in  Jeon’s  research  [10].  Table  (3.2) 
lists  the  r0  values  for  the  data  used  in  this  research. 

Depending  on  where  the  point  source  is  focused  on  the  CCD  array,  the  system  PSF 
will  change.  To  capture  these  changes  in  shape  of  the  system  PSF,  the  convolved  PSF 
must  be  shifted  and  then  resampled  at  increments  of  120  samples.  As  one  pixel  consists  of 
120  samples,  this  research  shifts  the  convolved  PSF  by  60  samples  up,  down,  left  and  right. 
This  is  consistent  with  the  Zingarelli,  et  al.  paper  and  table  (2.1).  However,  in  this  research, 
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Figure  3.4:  Zoomed  in  image  of  the  SST  non-downsampled  psf  (2400  x  2400  sample  size). 


Table  3.2:  Seeing  Parameters  for  the  Data  from  the  SST 


Day 

Seeing  parameter  (cm) 

73 

5.9 

74 

6.1 

75 

5.6 

81 

3.5 

82 

4.8 

83 

6.4 

the  shifts  are  done  in  a  different  order  which  is  reflected  in  Table  (3.3).  These  shifted 

system  PSFs  represent  the  nine  hypotheses  that  are  used  in  the  MHT;  these  hypotheses  are 
shown  in  Figure  (3.5).  Theoretically,  an  infinite  number  of  hypotheses  could  be  created  by 
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Table  3.3:  Alternative  Hypothesis  Sub-pixel  Shifts 


Hypothesis  (//,) 

Horizontal  Shift  (samples) 

Vertical  Shift  (samples) 

1 

-60 

-60 

2 

0 

-60 

3 

60 

-60 

4 

-60 

0 

5 

0 

0 

6 

60 

0 

7 

-60 

60 

8 

0 

60 

9 

60 

60 

dividing  the  pixel  into  smaller  and  smaller  shifts;  however,  in  this  research,  nine  hypotheses 
are  chosen  in  order  to  keep  computation  times  at  a  reasonable  level. 

3.3  Setup  of  Detection  Algorithms  in  Matlab 

The  two  main  detection  theories  used  in  this  research  are  the  point  detector  Binary 
Hypothesis  Test  (BHT)  and  the  Multi-hypothesis  Test  (MHT).  Both  these  detection 
algorithms  are  considered  matched- filter  tests  and  are  based  on  a  signal-to-noise-ratio 
(SNR)  that  is  compared  to  a  threshold  value,  r.  Essentially,  the  basic  equation  for  both 
tests  are  equivalent.  In  the  BHT,  it  is  a  ratio  of  whether  the  H0  hypothesis,  which  is  the 
likelihood  that  no  object  is  present  in  a  pixel,  as  compared  to  the  H\  hypothesis,  which  is  the 
likelihood  that  an  object  is  in  a  pixel.  The  MHT  takes  this  one  step  further  by  replacing  the 
H\  hypothesis  with  multiple  hypotheses,  each  of  which  are  compared  to  the  H0  hypothesis. 
Each  comparison  in  the  MHT  returns  a  SNR;  the  highest  valued  SNR  is  considered  the 
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HI 


10  20 
x-pixel 

H4 


10  20 
x-pixel 

H7 


10  20 
x-pixel 


10  20 
x-pixel 


H5 


10  20 
x-pixel 


H8 


10  20 
x-pixel 


H3 


10  20 
x-pixel 


10  20 
x-pixel 


10  20 
x-pixel 


Figure  3.5:  Nine  hypotheses  (H{  -  H9)  used  in  the  MHT  method 


most  likely.  As  each  hypothesis  in  the  MHT  is  created  by  shifting  the  focal  point  on  the 
CCD  array,  this  also  has  the  additional  benefit  of  providing  sub-pixel  position  information 
for  the  detected  object. 

3.3.1  Binary  Hypothesis  Test. 

The  baseline  point  detector  used  in  the  SST  can  be  derived  from  a  BHT;  it  is  expressed 
in  a  likelihood  ratio  test  (LRT)  as  follows: 

A=  P(d(x,yMx,y)e[l,Md]\Hi)>i 

P(d(x,yWx,y)e[l,Md]\H1)<  ’  ^  ’ 

Ho 

where  d(x,  y)  is  the  image  data,  x  and  y  are  the  pixel  coordinates  and  Md  is  the  number 
of  pixels  in  one  dimension  of  a  chosen  square  window  in  the  CCD  array  plane  [20].  The 
hypothesis  H\  is  a  likelihood  that  an  object  is  present  in  the  pixel  of  interest  and  H0  is  the 
hypothesis  that  an  object  is  not  present  in  the  pixel  of  interest.  P(d(x,  y)V(x,  y)  e  [  1 ,  Md\\Hi) 
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is  the  joint  conditional  probability  of  the  data  given  that  Hh  i  e  {0, 1}  is  true.  Therefore;  in 
the  BHT,  the  LRT  becomes  a  ratio  of  the  hypothesis  probabilities. 

Using  a  Gaussian  noise  distribution  to  represent  the  probability  P[d(x,y)\Hj]  given 
hypothesis  or  H(h  the  BHT  LRT  becomes 

FT  jnr  1  {-^[d(w^)-B-9h(w,z)]2)  R 

A  =  P[d{x,y)\Hx]  =  ,!=‘i  ^ _ > 

P[d(x,y)\H0\  <  ’ 

,  1  ,  dEncr 
W=  1  z=  1 

where  w  and  z  are  pixel  locations  in  the  window,  and  the  PSF  is  h(w,z )  [20].  The  joint 
conditional  probabilities  are  modeled  using  a  Gaussian  noise  distribution;  though  the  H\ 
joint  conditional  probability  would  be  more  accurately  modeled  by  a  Poisson  random  noise 
distribution  [20].  This  is  due  the  noise  which  limits  detection  performance  is  the  photon 
noise  of  the  stars,  optical  background,  and  target  [14].  As  the  amount  of  photons  arriving  at 
a  target  is  discrete,  yet  random,  it  is  most  accurately  modeled  using  a  Poisson  distribution. 
However,  as  research  in  [13]  showed  that  both  the  Gaussian  and  Poisson  algorithms  should 
have  close  to  the  same  level  of  performance  if  the  background  master  image  is  flat  and 
small;  this  allows  the  master  image  term  to  be  pulled  out  of  the  summation  and  distributed 
through  the  denominator.  This  results  in  the  same  expression  for  SNR  for  both  the  Gaussian 
and  Poisson  distributions  used  for  photon  statistics  [13].  As  the  SST  does  not  use  a  master 
background  image,  its  background  noise  distribution  can  be  considered  flat;  therefore,  a 
Gaussian  distribution  can  be  used.  The  value  Md  is  the  total  number  of  pixels  in  the  window 
of  the  image  [20].  The  value  B  is  the  background  noise  level,  which  is  computed  as  the 
local  sample  median  of  the  data.  The  6  is  the  space  object’s  irradiance  and  cr  is  the  standard 
deviation  of  the  noise  and  E\  ]  is  the  statistical  expectation  operation. 


B  =  median[d(w,z)\  V  (vv,  z)  e  [  1 ,  Md\ 


(3.19) 


cr  =  V E[d2(w,z )]  -  E[d(w,z)V 


\ 


Md  Md 

X  X  (d2(w,z)  -  B)2 

w=  1  z=  1 


(3.20) 
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To  simplify  Eq.  (3.18),  the  natural  logarithm  of  is  taken;  this  reduces  Eq.  (3.18)  to  the 
following 

Md  Md  Ih 

log{ Ag)  =  V  V  - [-2 B  ■  6h(w,z )  +  2 d(w,z)  ■  6h(w,z )  -  ( 9h(w,z ))2]  0.  (3.21) 

5„ 

This  can  further  rearranged  and  simplified  to 


Md  Md 

XX 


w=  1  z=l 


(J(w,  z)  -  B)h(w,  z)  >  6 

cr  <2 

Ho 


Md  Md 

XX 


W=  1  Z=1 


h2(w,z) 

cr 


(3.22) 


where  the  selection  of  d  will  be  chosen  to  achieve  the  desired  threshold  y  [20].  To  convert 
Eq.  (3.22)  into  an  equation  defined  in  terms  of  SNR,  a  new  random  variable  d2  must  be 
defined.  SNR  is  defined  as  the  ratio  of  the  mean  divided  by  the  standard  deviation.  The 
new  random  variable  d2  is  the  background  suppressed  data  and  must  have  a  zero  mean  in 
the  Ho  case.  The  random  variable  d2  is  defined  as 


d2(w,  z)  =  d{w,  z)  -  B. 


(3.23) 


The  correlation  of  the  system  PSF,  h(w,z),  with  the  background  suppressed  data  then 
becomes 


Md  Md 

XX  d2(w,  z)h(w  -  cx,  z  -  cy )  (3.24) 

W=  1  z=  1 

where  cx  and  cy  are  the  coordinates  of  the  pixel  being  tested  [20].  It  can  be  shown  that  the 
Ho  case  has  a  mean  /u2  of  zero  by 


Hi  =  E 


Md  Md 


XX  d2(w,  z)h(yv  -  cx,  z  -  cy) 


IW=\  Z=  1 
Md  Md 


(3.25) 


^  ^  E[d2(w,  z)\E[h(w  -  cx,  z  -  cy)]  =  0, 


W=  1  Z=1 
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and  a  variance,  cr2,  of 


a\  =  E 


=  E 


(  Md  Md  \21 

zz  Cx,  Z  -  Cy ) 

kW=l  Z=  1  t 

Mi  Mi  Mi  Mi 

zz  d2(w,z)h(w  -  cx,z  -  cy)  ZZ  d2(m,  ri)h(w  -  cx,  z  -  cy)  (3.26) 

w=\  z—  1  w=l  z—  1 

Mi  Mi  Mi  Mi 

=zzzz  £[d2(w,  z)]£[*/2(fra,  n)]/r(w  -  cx,  z  -  cy)h(m  -  cx,  n  -  cy ). 

w=  1  z—  1  /n=l  n=l 

Using  two  cases,  Eq.  (3.26)  can  be  solved.  The  first  case  is  when  w  ±  m  and  z  ±  n;  the 
second  case  involves  using  the  Dirac  delta  function,  6(w-m,z-n),  and  is  valid  when  w  =  m 
and  z  =  n.  Evaluating  Eq.  (3.26)  using  these  two  cases,  it  reduces  to 

Mi  Md 

(3.27) 

W=  1  Z=1 

The  standard  deviation  of  the  normalized  data  correlated  with  the  total  system  PSF  is  simply 
the  square  root  of  Eq.  (3.27).  Using  the  standard  deviation,  it  is  now  possible  to  express 
the  LRT  in  terms  of  SNR  [20] 

Mi  Md  Mj  Mi 


& 2  =  ^2YjYjh^W,Z)- 


SNRr 


Z  Z  (d(w,  z)  -  B)hcorr(w  -  cx,  z  -  Cy)  /-/,  f  Z  Z  (hcorr(w,  z))~ 

W—  1  1  >  VI’ =  i  .  I 


(T\I 

1  Mi  Mi 

Z  Z  h2corr(w,z) 

y 

w=  1  z=  1 

< 

Ho 


y.  (3.28) 


Mj  Md 


<J  A  Z  Z  h2corr(w,z) 


I  W=  1  Z- 1 

The  PSF  hcorr(w ,  z)  is  the  sampled  PSF  [20].  In  the  case  of  the  baseline  point  detector  BHT, 
hCorr( w, z)  is  one  pixel  represented  as  a  delta  function,  6(w  —  cx,z  —  cy);  this  reduces  Eq. 
(3.28) to 


Md  Mi 

Z  Z  (d(w,  z)  -  B)5(w  -  cx,  z  -  cy)  //, 

C..D  w=u=  i  \d(cx,  Cy)  -  B)  > 

^  N K baseline  ~  , -  —  T • 

Md  Mi  < 

Z  Z  d(w,z) 

V  w=  l  z=  l 


(3.29) 


3.3.2  Multi-Hypothesis  Test. 

There  are  two  forms  of  the  MHT  used  in  this  research.  Both  forms  rely  on  using  nine 
hypotheses,  Hi  -  Hg,  that  are  then  correlated  with  the  image  data.  The  only  difference 
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between  the  two  methods  is  in  the  calculation  of  the  standard  deviation  of  the  background 
noise.  In  the  basic  MHT,  the  noise  is  simply  the  standard  deviation,  cr,  of  the  background, 
B.  In  the  MHT  plus  method,  the  outliers  in  the  background  noise,  such  as  bright  stars,  are 
removed;  this  has  the  effect  of  lowering  the  background  noise  standard  deviation,  referred 
to  as  £  in  the  MHT  plus  method.  Being  in  the  denominator  of  the  MHT  equation,  this 
increases  the  SNR  in  cases  where  there  are  strong  noise  outliers  in  the  observed  image. 

3.3.2. 1  Basic  Multi-hypothesis  Test. 

The  SNR  equation  for  the  MHT  is 

Md  Md 

X  X (d(w, z)  -  B)hsamp(w  -  cx  —  aa,z  —  cy  —  J3a)  h„ 

SNRa  =  -  - >  y  MHT,  (3.30) 

Md  Md  < 

<rjz  E  h2  (w,  z) 

V  H’=l  Z=1 

where  aa  and  f3a  represent  the  half-pixel  horizontal  and  vertical  shifts  respectively  of  the 
sampled  PSF.  This  allows  the  MHT  to  provide  both  detection  information  as  well  as  sub¬ 
pixel  position  information.  The  MHT  is  in  contrast  to  the  BHT  which  only  tests  whether 
an  object  is  in  a  pixel  or  not.  Theoretically,  more  hypotheses  would  provide  a  greater 
level  of  sub-pixel  location  and  increase  detection  performance;  however,  a  finite  amount  of 
hypotheses  were  picked  in  order  to  make  the  test  numerically  practical.  In  performing  the 
MHT  method,  the  hypothesis  that  maximizes  the  SNR,  while  at  the  same  time,  increasing 
the  probability  of  detection  is  chosen.  The  choice  of  the  hypothesis  with  the  highest  SNR 
also  provides  sub-pixel  position  information  on  the  image  location  [20]. 

This  article  describes  three  detection  strategies  that  could  be  used  on  the  SST.  The  first 
technique,  which  is  currently  being  employed  on  the  SST,  is  the  baseline  point  detector 
[20].  The  second  detection  method  that  is  described  is  the  correlator  or  matched  filter 
detector.  Both  of  these  methods  are  considered  BHTs.  The  last  method  that  is  proposed 
is  the  MHT  method.  Of  these  various  methods,  the  MHT  is  the  method  that  is  proposed 
to  improve  the  detection  performance  of  the  SST.  Results  have  shown  that  the  probability 
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of  detection  by  using  MHT  improved  by  as  much  as  50  percent  over  existing  detection 
algorithms. 

The  two  BHT  methods  and  the  MHT  are  based  on  testing  one  hypothesis  against 
another  in  terms  of  SNR  [20].  The  first  hypothesis  is  known  as  H0  and  theorizes  that  a 
space  object’s  image  is  not  in  a  pixel.  The  alternative  hypothesis,  H\,  supposes  that  there  is 
a  space  object’s  image  in  a  pixel.  The  MHT  method  further  divides  a  pixel  into  sub-pixel 
hypotheses,  {H\  — //g),  and  compares  them  to  the  Ho  hypothesis.  The  MHT  method  gives  an 
advantage  by  mitigating  some  of  the  aliasing  caused  by  the  undersampled  SST  observation 
images.  In  the  MHT  method,  the  hypothesis  that  maximizes  SNR  while  increasing  the 
probability  of  detection  (Pd)  is  chosen  thereby  providing  sub-pixel  position  information 
while  increasing  the  Pd  over  both  BHT  methods  [20].  In  determining  the  precise  location  of 
a  space  object,  from  which  angular  rate  can  be  determined  over  subsequent  image  frames, 
the  sub-pixel  location  provided  by  the  MHT  method  provides  a  significant  advantage  over 
the  other  detection  methods. 

The  MHT  is  favored  over  the  correlator  detector  because  the  shape  of  the  sampled  PSF 
changes  depending  on  where  the  object  is  imaged  on  the  observation  pixel  array  [20];  the 
correlator  detector  does  not  take  this  into  account.  In  the  MHT,  the  sampled  PSF  changes 
depending  on  the  location  in  the  pixel  array.  These  various  sampled  PSFs  give  rise  to  the 
various  hypotheses  in  the  MHT  detector.  The  SNR  for  the  MHT  is  defined  as 

c  . TD  Ijwd  Hfrf(d(w, z)  -  B)hsamp(w  +  cx-aa,z  +  cy- (3a)  > 

SNRa  -  -  - : -  jM-ary  (3.31) 

CT  '  >/ Xwd  Hf d  h2samp{w,  Z )  Ho 

In  equation  (3.31),  aa  and  [3a  are  the  horizontal  shift  and  vertical  shift  respectively  that 
correspond  to  a  sub-pixel  location;  the  nine  hypotheses  are  derived  from  combinations  of 
+  15 /urn  and  0 tun  shifts  in  aa  and  [3a  [20].  These  hypotheses,  hsamp,  are  then  correlated 
with  the  image  data,  d{w,  z)  minus  the  background  noise,  B.  The  square  root  portion  of 
the  denominator  in  equation  (3.31)  normalizes  the  correlation.  The  cr  denotes  the  noise 
standard  deviation.  Each  SNR„  corresponds  to  a  hypothesis,  Ha\  from  this,  it  is  deduced 
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that  the  maximum  SNR„  corresponds  to  the  hypothesis  that  most  accurately  reflects  the 
sub-pixel  location.  The  maximum  SNR„  is  considered  the  true  SNR  of  the  detected  object. 

S  NR-M-ary  =  max(S  NRa)  (3.32) 

The  ~y M-ary  denotes  the  threshold  value  that  the  SNRW_,„.V  must  be  greater  than  in  order  for 
an  object  to  be  considered  as  detected  [20].  Mathematically,  a  low  threshold  value  increases 
the  probability  of  false  alarm,  PFA;  for  this  reason,  the  threshold  value  is  generally  increased 
in  order  to  minimize  the  PFA.  For  SST,  the  threshold  value  is  set  to  6.2212  so  that  the  PFA 
is  9.87e-10. 

By  using  the  MHT,  it  is  hoped  that  both  the  detection  of  earth  orbital  objects  will 
increase  as  well  as  provide  additional  sub-pixel  location  information  by  which  more 
accurate  angular  rates  of  orbital  objects  can  be  made. 

3.3.2. 2  Multi-hypothesis  Test  with  outlier  removal. 

The  Zingarelli,  et  al.  paper  also  proposes  a  separate  method  for  computing  the 
standard  deviation,  cr,  of  the  background  noise  as  listed  in  Eq.  (3.20).  This  method  rejects 
any  noise  sample  in  the  window  surrounding  the  pixel  to  be  tested  whose  values  are  outside 
those  predicted  by  Gaussian  statistics  [20] .  This  has  the  effect  of  eliminating  extremely  high 
noise  or  nearby  stars  in  the  calculation  of  background  noise.  In  this  alternative  method,  the 
background  noise,  B ,  is  still  computed  as  in  Eq.  (3.19).  The  squared  deviations,  D,  from 
the  background  within  the  window  are  computed  by 

Dim)  =  (d{m)  -  B)2,  (3.33) 

where  m  is  a  pixel  in  the  window  around  the  pixel  being  tested  [20].  The  mean,  M,  of  the 
squared  deviation  D  is 


M  =  E[D(m)\, 


(3.34) 
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with  a  standard  deviation  of  S  computed  as 


S  =  V E[D2(m )]  -  E\D(m)\2 


Md 

X  D2(m) 


A 


m=  1 

“^5 


-M2. 


(3.35) 


The  new  noise  standard  deviation,  £,  is  then  computed  by  excluding  any  pixel,  m,  that  is 
outside  three  standard  deviations  from  the  mean;  mathematically,  this  excludes  any  pixel,  m 
where  Dim)  >  ( M  +  3-S ).  The  new  standard  deviation  implemented  in  the  MHT  has  shown 
to  improve  detection  performance  by  lowering  the  normalizing  factor,  thus  increasing  the 
SNR.  This  is  implemented  in  Eq.  (3.30)  by  replacing  cr  with  £  [20]. 


3.4  Setup  of  Position  Algorithms  in  Matlab 

The  position  algorithms  used  in  this  research  are  based  on  the  results  of  the  BHT  and 
MHT  detection  algorithms.  There  were  three  types  on  position  algorithms  used  in  this 
research.  The  first  method,  referred  to  as  the  centroid  method,  takes  the  results  of  the  BHT, 
which  are  ones  and  zeros,  and  centroids  these.  In  most  cases,  the  BHT  identifies  a  single 
pixel  as  a  detection.  In  this  case,  the  integer  value  of  the  pixel  is  used  as  the  position  of  the 
detected  object.  This  position  algorithm  is  currently  the  method  employed  on  the  SST.  The 
next  method  used  in  this  research  is  the  Center  of  Intensity  (Col)  algorithm.  This  method 
maps  the  detected  points  back  to  the  original  data,  d(cx,cy).  The  intensity  values  around 
and  including  the  detected  point  are  then  centroided  to  get  a  more  accurate  position.  The 
last  method  uses  the  SNR  values  from  the  MHT.  In  this  method,  four  of  the  hypotheses’ 
SNR  values  are  recorded  in  a  master  array.  A  centroid  operation  is  then  performed  on 
the  master  array  to  determine  the  position.  This  method  has  calculated  the  most  accurate 
position  information  for  a  detected  space  object. 

3.4.1  Centroid  Algorithm. 

The  detection  of  an  object  in  the  BHT  is  done  by  simply  noting  in  an  array  whether  an 
object  is  detected  (set  to  a  value  of  1)  or  not  detected  (set  to  a  value  of  0)  in  the  detection 
array,  D(x,y).  Occasionally,  a  detected  object  is  bigger  than  a  single  pixel.  In  this  case,  a 
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2  4  6  8  10  12  14  16  18  20 

x-pixel 

Figure  3.6:  Image  of  star  viewed  by  the  SST. 


centroid  algorithm  is  used  to  determine  its  pixel  location.  The  centroid  algorithm  for  xcenter 
(y center  is  identical)  coordinate  in  the  detection  array,  D(x,  y),  of  20  by  20  pixels: 

I"  0(*oO-* 

W=  (  } 

The  detection  array  is  simply  an  array  of  ones  and  zeros  marking  the  detection  of  objects 
with  the  value  one.  The  location  in  the  array  of  the  value  corresponds  to  the  pixel  location 
in  the  actual  image  frame.  Figure  (3.6)  shows  a  star  observed  by  the  SST.  The  BHT  run 

on  this  particular  frame  of  data,  shows  that  four  center  pixels  each  met  the  point  detection 
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Figure  3.7:  Detection  array,  D(x,y),  from  the  BHT  of  star  viewed  by  the  SST. 


threshold  and  thus  were  marked  with  a  value  of  one;  the  results  of  the  BHT  are  shown 
in  Figure  (3.7).  Using  Eq.  (3.36)  on  this  figure  would  yield  a  position  of  (11.5, 11.5)  as 
the  position  of  this  detected  space  object.  Looking  at  Figure  (3.6),  the  brightest  pixel  is 

(12. 12) .  It  is  assumed  in  this  research  that  the  actual  position  of  the  object  is  closer  to  the 

(12. 12)  pixel;  however,  as  the  point  detector  BHT  assigns  an  equal  weighting  to  detection 
points  and  does  not  take  into  account  the  intensity  values,  that  there  is  an  inherent  error  in 
the  position  measurements  of  the  centroid  method. 


39 


Table  3.4:  4x4  array  of  intensity  values,  dssr(x,y).  In  this  example,  the  position  of  the 
intensity  values  would  would  align  with  the  detection  points  in  Figure  (3.7). 
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3.4.2  Center  of  Intensity  Algorithm. 

The  center  of  intensity  takes  the  centroid  operation  one  step  further.  Instead  of 
performing  the  centroid  operation  on  the  detection  array 

D(x,y) 


,  the  centroid  operation  is  performed  on  the  actual  intensity  values  of  the  SST  data, 
dssr(x,y).  The  array  of  intensity  values  is  determined  by  multiplying  detection  array  by 
the  original  data  frame  from  the  SST;  this  effectively,  multiplies  by  zero  anything  that  does 
not  meet  the  threshold  of  the  BHT.  Figure  (3.8)  shows  an  image  of  dSST{x,y)  after  it  has 
been  multiplied  by  the  detection  array  D(x,  y).  The  centroid  operation  is  performed  using 
Eq.  3.37.  Table  (3.4)  contains  the  intensity  values  of  the  dssr(x,y )  that  are  used  in  the 
centroid  example  in  section  3.4.1. 

ZLZv=i  dSST(x,y)-x 

x center  -  3  3  \5.5t) 

Xjc=i  Tty=i  dssr(x,y ) 


Using  the  center  of  mass  equation  for  the  data  in  table  (3.4),  the  center  of  intensity  is 
(1 1.641, 11.567).  Based  on  an  observation  of  the  intensity  pattern  figure  (3.8),  it  should  be 
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2  4  6  8  10  12  14  16  18  20 
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Figure  3.8:  The  BHT  detection  points  mapped  back  to  intensity  values  of  star  viewed  by 
the  SST. 

more  clear  that  the  position  of  the  detected  object  is  weighted  and  closer  to  the  brightest 
pixel,  which  at  (12, 12). 

3.4.3  Center  of  SNR  Values  Algorithm. 

The  center  of  SNR  (CoS)  relies  strictly  on  the  results  of  the  MHT  method.  Each 
hypothesis  is  correlated  with  the  underlying  SST  data.  In  the  correlation  operation,  the 
PSF  is  slid  through  data  and  each  pixel  returns  a  SNR  value.  Some  of  the  hypotheses 
in  the  MHT  have  identical  intensity  patterns,  but  are  shifted  in  position.  For  example,  in 
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Figure  3.9:  The  four  hypotheses  of  the  MHT  mapped  to  the  SNR  master  matrix. 


Figure  (3.5),  H\,H3,H1 ,  and  H 9  are  identical  in  their  intensity  patterns,  but  are  shifted 
in  position.  This  is  true  for  H4  and  H 6,  as  well  as,  H2  and  H 8.  Due  to  the  correlation 
operation  sliding  the  PSF  through  the  data,  these  identical  PSF  intensity  patterns  will 
produce  identical  correlation  values  with  the  SST  data.  For  this  reason,  only  H 5,  H 6, 
H 8,  and  H 9  are  used  in  this  research. 

To  determine  the  position  using  the  CoS  method,  the  four  hypotheses’  SNR  values  are 
mapped  to  a  matrix  that  is  double  the  size  of  the  original  data  window.  Figure  (3.9)  shows 
the  SNR  master  matrix  values  for  same  example  used  in  the  centroid  and  Col  sections. 

The  position  of  the  object  is  found  by  performing  a  centroid  operation  on  the  SNR  master 
matrix  and  dividing  by  two,  which  accounts  for  the  SNR  master  matrix  being  twice  the 
size  of  the  original  data  matrix.  Performing  the  centroid  operation  on  this  example  would 
give  a  position  of  (3.2554,3.1771)  on  the  6x6  SNR  matrix  and  (1.6277, 1.5886)  on  the 
original  3  x  3  set  of  observation  data.  Each  of  the  four  hypotheses  correspond  to  a  sub-pixel 
position.  The  sub-pixel  position  of  the  object  is  then  weighted  by  each  of  the  SNR  values 
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Position  (x) 


Figure  3.10:  The  difference  in  position  between  the  MHT  position  algorithms  and  the 
centroid  and  Col  algorithms. 


from  each  of  the  hypotheses.  In  this  example,  the  H 9  hypothesis  has  the  strongest  SNR 
values.  The  H9  hypothesis  corresponds  to  the  bottom,  right  corner  of  a  pixel.  Therefore,  in 
this  example,  the  sub-pixel  position  is  weighted  most  strongly  to  the  H 9  sub-pixel  position. 
This  approach  is  performed  with  both  the  MHT  and  MHT  plus  methods. 

It  is  worth  noting  that  the  MHT  methods  introduce  a  bias  as  a  result  of  combining  all 
the  SNRs  into  one  master  matrix.  Figure  (3.10)  shows  an  example  of  this  bias  in  position 
over  several  frames  of  data  between  the  MHT  position  algorithms  and  the  centroid  and  Col 
position  algorithms.  As  the  angular  rate  is  measured  by  the  change  in  pixel  position  from 
frame  to  frame,  the  actual  position  information  does  not  affect  the  angular  rate  calculations; 
however,  in  this  research,  a  constant  is  added  to  the  MHT  position  results  to  bring  the 
position  back  in  line  with  the  Col  method. 

3.5  Measuring  Angular  Rate  from  Position  Data 

Measuring  the  angular  rate  from  the  position  data  requires  at  least  two  sequential 
frames  of  data  where  a  space  object  is  detected.  Once  the  position  is  known  in  the  two 
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frames  of  data,  then  it  becomes  a  mathematical  operation  to  measure  the  change  in  position 
over  a  period  of  time.  As  the  range  of  a  detected  space  object  is  generally  unknown,  only  the 
angular  rate  of  the  detected  object  can  be  surmised.  Each  pixel  in  the  focal  array  represents 
an  angle  in  radians  of  the  field-of-view. 

3.5.1  Determining  Velocity  from  Position  Data. 

The  next  step  in  this  project  after  detecting  an  object  is  to  determine  an  object’s 
velocity  in  pixels  per  second.  To  do  this,  it  involves  measuring  the  change  in  position 
(pixel  location)  between  two  subsequent  frames  of  the  detected  object.  In  this  research,  the 
changes  is  in  the  x  pixel  position,  represented  as  AY  and  the  changes  in  y,  denoted  as  Av, 
are  done  separately.  Figure  (3.11)  shows  an  example  of  the  calculated  angular  rate  in  the 
x  -  coordinate  direction;  note  that  the  angular  rate  in  this  figure  is  in  radians  per  frame  of 
data.  This  is  a  simple  measurement  of  change  between  two  positions  as  described  by: 


A.v  =  x2  -  xu 


(3.38) 


Ay  =  >’2  -  >’l  • 

From  the  total  magnitude  of  the  change  in  pixel  position,  Axv  is  also  determined  using  the 
following  equation: 


A  xy  ~ 


(3.39) 


In  this  research,  only  the  magnitude  Axy  is  considered.  While  the  direction  portion  of 
the  change  vector  is  important  for  determining  the  correct  orbit,  this  research  is  focused 
on  determining  which  angular  rate  algorithm  is  most  accurate.  The  measurement  of  the 
angular  rate  of  a  star,  which  has  a  well  known  and  precise  is  used  to  determine  which 
angular  rate  algorithm  is  most  accurate. 


3.5.2  Transforming  Change  in  Pixel  Position  into  Angular  Rate. 

In  order  to  determine  the  angular  rate,  it  is  necessary  to  determine  what  each  pixel 
represents  in  terms  of  radians.  This  is  done  by  taking  the  tangent  of  the  size  of  each  pixel, 
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radians/frame  divided  by  sidereal  rate  (x) 


Figure  3.11:  Example  of  calculated  angular  rate  in  radians  per  frame 


psize,  in  the  SST  divided  by  the  focal  length,  /,  of  the  SST.  The  pixel  size  of  the  SST 
is  15 pm  x  1 5pm;  however,  the  pixels  are  binned  together  in  a  two  by  two  set  effectively 
creating  a  30 pm  x  30 pm  pixel.  The  focal  length  of  the  SST  is  3.5 m.  Thus  the  equation  for 
determining  what  each  pixel  represents  in  radians  for  the  SST  is: 

Ae  =  tan|^y^j  =  tan  | j  -  8.571  prad/ pixel.  (3.40) 

The  simplest  expression  for  determining  the  angular  rate  is  Ae  divided  by  the  interval  of 
time,  t,  over  which  the  change  occurs.  The  frames  of  data  used  in  this  research  from  the 
SST  are  taken  at  1  second  intervals.  The  angular  rate,  9,  of  an  objects  is  then: 

Ae-  A xy  radians 

9  -  -  (5.41) 

t  second 

Figure  (3. 1 1)  is  an  example  of  the  total  magnitude  of  the  angular  rate  of  a  star  measured  in 
radians  per  frame. 
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3.6  Determination  of  Statistics 

The  method  for  comparing  which  angular  rate  algorithm  is  the  most  accurate  in 
determining  the  angular  rate  is  based  on  statistical  analysis.  For  this  research,  stars,  which 
have  a  highly  precise  and  constant  angular  rate,  are  used  as  a  way  to  compare  a  calculated 
angular  rate  with  a  known  angular  rate. 

Stars  appear  to  make  one  revolution,  2 n  radians,  of  the  Earth  per  sidereal  day  which 
is  23.9344699  hours  [6].  To  calculate  the  angular  rate  in  radians  per  second,  denoted  as  dT, 
where  T  represents  the  true  angular  rate.  The  true  angular  rate  of  a  given  star  is 

BT  -  — — — -  - =  72.921 urad/sec.  (3.42) 

(23.9344699/ms  •  3600.vec) 

To  simplify  the  interpretation  of  the  data,  the  calculated  angular  rate  was  divided  by  the 
true  angular  rate  of  a  star.  As  only  stars  are  used  in  this  research,  the  ratio  of  calculated 
angular  rate  divided  by  the  true  angular  rate  should  be  approximately  equal  to  one.  Figure 
(3.12)  shows  an  example  plot  of  this  ratio;  as  can  be  seen  in  the  plot,  the  angular  rate  hovers 
around  1. 

There  are  two  ways  to  compare  the  three  angular  rate  algorithms,  centroid,  Col,  and 
CoS,  used  in  this  research.  The  first  method  involves  determining  the  average  of  the  angular 
rate  measurements.  Mathematically,  this  is  simply  the  expected  value,  £[],  of  all  calculated 
angular  rates,  E\X\,  where  X  is  the  data  set  consisting  of  the  angular  rate  measurements  at 
each  point,  (6N-\,6N-2, #2,  #1),  and  N  represents  the  number  of  sequential  frames  where 
a  star  is  detected.  The  angular  rate  algorithm  with  expected  value  closest  to  a  value  of 
one  can  then  be  considered  the  most  accurate  method.  A  concern  about  using  this  method 
is  that  there  can  be  large  variances  in  the  calculated  angular  rate  while  still  having  an 
expected  value  close  to  one.  The  SST  builds  its  tracks  for  a  detected  object  from  only  three 
detections;  therefore,  the  variance  of  the  angular  rate  becomes  a  significant  factor. 
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Figure  3.12:  Example  of  calculated  total  magnitude  of  angular  rate  in  radians  per  frame 
divided  by  sidereal  rate 

To  calculate  the  variance  of  the  angular  rate,  the  standard  mathematical  expression  for 
variance  is  used;  the  variance  is 


(3.43) 


The  angular  rate  algorithm  that  has  the  smallest  cr2  can  be  considered  the  most  consistent 
algorithm. 

The  comparison  of  angular  rate  algorithms  can  then  be  made  on  two  bases.  First,  the 
output  of  the  most  accurate  angular  rate  algorithm  can  be  made  on  the  comparison  with 
the  true  angular  rate,  0T,  by  determining  the  expected  value  of  the  observation  data  set, 
X.  Second,  a  judgment  on  the  consistency  of  measurement  of  the  angular  rate  calculation 
based  on  the  variance  of  the  observation  data  set  can  be  measured. 

3.7  Error  Estimation 

The  two  primary  sources  of  phase  errors  on  a  light  wave  as  it  propagates  from  the 
source  in  space  to  the  SST  are  atmospheric  turbulence  and  the  optics  of  the  SST  [15]. 
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These  phase  errors  in  the  wavefront  have  an  effect  on  measuring  the  position  of  a  detected 
object.  As  the  Zernike  coefficients  used  to  model  the  SST  are  based  on  a  mean  of  36 
DONUT  algorithms  run,  it  is  assumed  that  these  are  most  accurate  method  to  determine  the 
Zernike  coefficients,  therefore  this  section  will  not  discuss  this  type  of  phase  error  in  more 
detail.  The  tilt  phase  aberration  due  to  atmospheric  turbulence,  which  causes  an  object 
to  be  shifted  in  position  in  the  observation  plane,  is  the  most  significant  error  associated 
with  atmospheric  turbulence;  therefore,  this  section  will  deal  primarily  with  analyzing  the 
effects  due  to  tilt.  The  error  associated  with  making  measurements  of  position  using  an 
optics  system  is  called  the  observation  uncertainty.  The  first  part  of  this  section  will  deal 
with  determining  the  observation  uncertainty  for  the  SST  based  on  various  tilt  phase  errors. 
The  second  part  of  this  section  will  deal  with  how  the  observation  uncertainty  affects  the 
measurements  of  angular  position,  rate,  and  acceleration. 

3.7.1  Lower  Bound  on  Position  Error  Due  to  Atmospheric  Tilt. 

Atmospheric  tilt  phase  error  affects  the  performance  of  the  SST  by  changing  the 
geometric  position  of  the  target  on  the  detector  plane.  The  long  exposure  tilt  is  often  treated 
as  a  Gaussian  random  variable  with  a  zero  mean  and  unit  variance  [15].  Figure  3.13  shows 
an  example  of  the  Gaussian  shape  and  blurring  effects  on  a  long  exposure  image.  However, 
at  any  given  instant  in  time,  tilt  only  causes  a  translation  or  shift  in  the  geometric  position 
of  the  observed  object.  The  observed  images  at  these  instants  in  time  are  often  referred 
to  as  short  exposure  images.  Figure  3.14  shows  an  example  of  the  diffraction  limited  PSF 
generated  by  the  SST.  Notice  that  the  PSF  is  slightly  off  the  center  point  of  (600,  600)  in 
the  image;  it  is  actually  closer  to  the  (609,  600).  The  short  exposure  PSF  example  shows  a 
translation  due  to  atmospheric  turbulence,  but  no  blurring  due  to  atmospheric  turbulence. 

These  data  used  in  this  project  from  the  SST  have  an  exposure  or  integration  time  of 
100ms.  This  integration  time  has  characteristics  of  both  a  long  and  short  exposure  image. 
In  order  to  model  SST,  the  100  ms  integration  time  is  divided  evenly  into  1  ms  increments. 
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Figure  3.13:  Long  exposure  with  r0  =  15cm 


Figure  3.14:  Short  exposure  with  r0  =  15  cm 


A  phase  screen  is  generated  for  each  1  ms  increment  and  a  corresponding  PSF  is  generated 
from  each  phase  screen.  The  1  ms  increment  is  chosen  in  order  to  eliminate  the  effects  of 
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any  time  averaging  of  the  image  due  to  the  atmosphere  continually  and  rapidly  changing 
with  time.  The  general  guideline  that  a  short  exposure  image  is  somewhere  between  0.01 
and  0.001  seconds  or  shorter  depending  on  the  wind  velocity  [8]. 

Fried’s  seeing  parameter,  r0,  is  used  to  generate  each  phase  screen  and  was  considered 
to  stay  constant  over  the  100  ms  integration  time.  Each  phase  screen  that  is  generated  is 
statistically  correlated  with  the  preceding  phase  screen  in  order  to  accurately  model  the 
temporal  behavior  of  the  atmosphere  [15].  The  atmospheric  phase  screen  is  assumed 
to  evolve  in  time  via  Taylor’s  frozen  flow  hypothesis  [15].  Taylor  hypothesized  that 
temporal  variations  of  turbulent  eddies  at  a  point  are  produced  by  the  horizontal  transport 
of  atmospheric  properties  of  these  turbulent  eddies  by  the  mean  wind  speed  flow  and  not 
by  changes  in  the  turbulent  eddies  themselves  [5].  The  time  scale  associated  with  temporal 
variation  of  the  turbulent  eddies  due  to  wind  is  typically  on  the  order  of  1  second,  whereas 
the  temporal  variations  of  the  eddy  turnover  is  much  slower,  typically  on  the  order  of  10 
seconds  [5].  The  effect  of  the  wind  on  the  turbulent  eddies  is  to  translate  the  phase  screen 
a  distance  in  the  horizontal  and  vertical  directions  equal  to  the  wind  velocities,  vx  and  vy, 
multiplied  by  the  time  between  the  sample  images  At  [15].  In  this  research.  At  was  1  ms. 
The  relationship  between  the  phase  errors  and  wind  speed  at  times  0  and  t2  under  Taylor’s 
frozen  flow  hypothesis  is  equal  to  Eq.  (3.44)  [15]  where  the  atmosphere  is  considered 
statistically  homogeneous  [5]. 


6atm(w,  S,  h  +  At)  =  0atm(w  +  VXA t,  S  +  VyAt,  tj)  (3.44) 

The  statistical  correlation  between  each  phase  screen  is  defined  by  the  phase  structure 
function  for  Kolmogorov  turbulence 

,T2x  +  T2  5/6 

De(r,,r,)  =  6.88(-^)  (3.45) 

rQ 
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where  rr  =  w2-wq  and  rv  =  s2-s\  represent  the  spatial  frequencies  in  the  x  and  y  directions 
respectively.  Equation  3.45  is  then  used  to  determine  the  phase  structure  function 


De(At) 


Ty)[De(TX  +  VXA t,  Ty  +  VyAt)  ~  D 0(T x,  Ty)\dTXdTy 


(3.46) 


where  P,  is  the  aperture  function  [15].  The  phase  structure  function  De(Al)  is  then  used  to 
compute  the  tilt  correlation  DA  At  )  through  the  following  relationship: 


De(At)  =  2fle(0)  -  2Re(At). 


(3.47) 


RAO)  is  equal  to  the  variance  of  the  horizontal  tilt  07  and  the  vertical  tilt  crj,.  Both  07  and  cr2 
are  considered  equal  and  are  related  to  the  aperture  diameter  Dr  and  the  seeing  parameter 
r0  through  the  following  equation  [15]: 

oo  /  Dr  \ 5/3 

cr2e  =  o-l  =  RAO)  =  0A4S(—)  .  (3.48) 

v  r0  ' 


Equation  3.48  is  the  Zernike-Kolmogorov  residual  error  for  the  tilt  phase  errors  e  and  /J 
developed  by  Noll  [12].  RAAt )  is  the  temporal  autocorrelation  and  is  the  expected  value 
of  the  product  of  the  tilt  at  two  different  times.  This  is  mathematically  expressed  as: 
RAAt)  =  E[e{t\)e{t2)\  [15].  Once  RAAt)  is  known,  the  conditional  probability  of  the  tilt  at 
time  t\  +  At  given  the  tilt  at  time  t{  is  a  Gaussian  random  variable  with  a  mean  equal  to  [15] 


£[e2|ei]  =  ei 


RAAt) 


cr- 


(3.49) 


The  conditional  variance  is  for  the  Gaussian  random  variable  is 


RAAt)  A 


CT^ 


)  I6' 


crt-Rl(At) 


cr- 


(3.50) 


Once  the  conditional  mean  and  variance  are  known,  new  phase  screens  can  be  generated. 
The  new  conditional  tilt  value  e2  (the  same  equations  apply  to  J32)  is  given  by 


RAAt) 

e 2  =  61 - +  creg2 

Cl 


(3.51) 
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where  g2  is  a  normally  distributed  pseudorandom  number  [15].  The  new  tilt  phase  error  is 


0atm2(u,  s)  =  e2(u/Dr)  +  J32(v/Dr).  (3.52) 

The  atmospheric  transfer  function  Hatm(u,v)  is  then  equal  to  e'e“,m^<u'v) .  The  total  transfer 
function  is  Htot(u,v )  =  Hatm(u,v)Hoptics(u,v)  [15].  Hoptics(u,v)  is  the  optical  transfer 
function  of  the  SST  and  includes  the  SST’s  optical  aberrations  defined  by  Zemike 
polynomials  2  through  12.  Htot(u,v )  is  then  used  to  create  a  new  random  conditional  PSF 
for  each  1  ms  increment. 

In  this  research,  only  one  phase  screen  per  1  ms  increment  is  used  to  model  the  effects 
of  the  atmospheric  turbulence.  Propagation  along  a  slant  or  vertical  path  requires  an  index 
of  refraction  structure  parameter,  C2(h),  that  changes  as  a  function  of  altitude  (h)  [5].  The 
size  of  an  object  and  the  distance  that  the  light  propagates  from  that  object  forms  an  angle 
6.  If  this  angle  is  less  than  the  isoplanatic  angle,  then  the  atmospheric  turbulence  can  be 
modeled  with  a  single  phase  screen  because  the  angle  over  which  atmospheric  turbulence 
is  modeled  is  essentially  unchanged  [5].  One  of  the  most  popularly  used  atmospheric 
turbulence  models  to  express  the  C2(h)  parameter  is  the  Hufnagel- Valley  (H-V)  model. 
This  is  mathematically  described  in  equation  3.53  where  h  is  in  meters  [m],  w  is  the  rms 
windspeed  in  m/ s  and  A  is  a  nominal  value  of  C2(0)  at  the  ground  in  m~2/ 3  [5]. 

C2n(h)  =  0.00594(w/27)2(10“5/r)1  V/!/10()0  +  2.7  •  10"16e"*/150°  +  Ae~hlwo  (3.53) 

An  often  used  version  of  Hufnagel- Valley  is  determined  by  setting  w  =  21  m/s  and 
A  =  1.7  •  10 ~l4m~2/3;  this  is  commonly  referred  to  as  H-V5/7  [5].  For  a  wavelength 
A  =  0.5/rm,  the  5  in  H-V5/7  means  that  a  value  of  5cm  is  calculated  for  the  atmospheric 
coherence  diameter  and  the  7  refers  to  a  7  /rrad  isoplanatic  angle  [5].  Considering  H-V5/7, 
about  the  largest  object  that  can  be  observed  at  geosynchronous  orbit  (36,000  km)  and  still 
meet  the  criteria  under  H-V5/7  for  the  7  /rrad  isoplanatic  angle  is  252.0m;  see  Figure  3.15. 
It  can  be  assumed  that  most  of  the  objects  observed  by  the  SST  are  smaller  than  252.0m 
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and  therefore  meet  the  criteria  of  being  isoplanatic.  Therefore,  the  use  of  only  one  phase 
screen  to  model  the  phase  error  from  atmospheric  turbulence  is  valid. 


36,000,000  m 


t 


h252.0m 


Figure  3.15:  Isoplanatic  Angle,  6.  Note,  the  figure  is  not  to  scale. 


The  final  step  to  create  a  representative  100ms  exposure  image  from  SST  is  to  average 
all  100  short  exposure  1  ms  PSFs  together  to  create  a  mean  PSF.  The  mean  PSF  is 
considered  to  be  an  accurate  representation  of  an  image  generated  by  the  SST  over  the  100 
ms  integration  time.  The  mean  PSF  is  then  downsampled  and  integrated  using  Eq.  (3.6) 
over  each  pixel  to  create  a  sampled  PSF.  Figure  3.16  shows  an  example  of  a  simulated  100 
ms  image  by  SST. 

In  order  to  measure  the  amount  of  translation  for  each  mean  PSF,  the  center  of 
intensity  algorithm  described  in  Section  3.4.2  is  implemented  to  find  the  position  of  the 
sampled  PSF.  The  translated  position  is  then  compared  to  a  ’’true”  PSF  where  no  tilt 
was  introduced  and  the  magnitude  of  its  difference  is  recorded.  To  find  the  magnitude 
of  the  difference,  G,  between  the  two  center  of  intensity  points,  a  simple  vector  magnitude 
equation  was  used  as  described  in  Eq.  (3.54). 


(3.54) 


The  magnitude  of  the  difference  could  then  be  thought  of  as  an  the  error  in  observed  pixel 
location  from  the  true  pixel  location. 
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The  phase  screen  generated  for  each  r0  is  generated  by  a  statistical  random  process. 
It  is  possible  to  find  a  statistical  mean  and  variance  for  this  random  process.  To  do  the 
statistical  analysis,  30  repetitions  are  run  for  a  particular  r0  and  wind  velocity  and  the 
magnitude  of  position  error  is  recorded  for  each  repetition.  The  mean  value  of  the  position 
error  is  expected  to  be  near  0  for  each  r0  value  as  the  tilt  phase  error  is  random  process  that 
is  based  on  a  zero  mean  Gaussian  distribution.  The  standard  deviation  of  the  position  error 
over  the  30  repetitions  is  the  key  statistical  parameter  that  is  trying  to  be  determined  in  this 
research.  The  standard  deviation  can  be  thought  of  as  an  error  in  pixel  position  for  a  given 
r0  and  wind  velocity  for  the  SST. 

In  this  research,  394  simulations  were  run  without  Poisson  noise.  Another  197 
were  run  with  the  inclusion  of  Poisson  noise  added  to  the  sampled  PSF.  The  numbers 
of  simulations  run  was  restricted  to  the  amount  of  time  available  for  research.  In  both 
simulations,  the  equation  representing  the  best  line  fit  did  not  change  after  approximately 
100  simulations.  Additionally,  The  inclusion  of  Poisson  noise  represents  a  more  accurate 
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model  as  the  arrival  of  photons  is  a  random  event;  therefore,  on  average  there  is  an  expected 
number  of  photons  arriving  at  any  given  time  interval.  However,  the  actual  amount  arriving 
at  any  given  moment  is  best  described  by  the  Poisson  impulse  process  [8].  If  K  is  the 
random  variable  representing  the  number  of  photons  arriving  at  a  given  time  interval  and 
A(t)  represents  the  rate  function  or  mean  number  of  photons  arriving  for  a  time  interval, 
then  the  probability,  P(K;  / , ,  t2),  that  K  photons  fall  within  the  time  interval  (t\  <  t  <  to)  is 
given  by  [8] 


P(K;  t\, 


\K 


VI 


f  A(t)dt  exp  f  A(t)dt 


VI 


K\ 


(3.55) 


Each  simulation  performed  30  repetitions  for  a  specific  r0,  wind  x-direction,  and  wind  y- 
direction  values  generated  from  a  Gaussian  random  number  generator  in  order  to  provide 
statistical  significance. 

The  results  of  these  simulations  showed  an  decreasing  error  in  position  as  r0  values 
increased;  these  results  are  shown  in  figure  (3.17)  and  (3.18).  In  figure  (3.17)  and  (3.18), 
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♦  magstdev 
—  Power  (magstdev) 


Figure  3.17:  Position  error  ±  I  d)error  vs.  r0 


the  y-axis  is  ±Oerror  in  units  of  pixels;  the  x-axis  is  the  rQ  values  in  units  of  centimeters.  A 
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♦  magstdev 
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Seeing  pararmeter,  r0(cm) 


Figure  3.18:  Position  error  +  \<frerror  vs.  r0 


best  line  fit  in  Microsoft  Excel  was  used  on  the  ±  1  standard  deviations  line  (marked  with 
diamonds  in  figure  3.17).  The  best  line  fit  equation  is: 

+  10error  =  0.0542r-°-827  (3.56) 

and 

±10error  =  0.071  lr-°-742  (3.57) 

for  the  model  with  Poisson  noise.  It  is  interesting  to  note  that  the  inclusion  of  Poisson 
noise  actually  reduced  the  error  due  to  the  Poisson  function  eliminating  small  values,  thus 
reducing  the  size  of  the  sampled  PSF.  Eqs.  (3.56)  and  (3.57)  states  that  for  a  particular  r0 
value  that  the  position  of  the  object  has  an  error  of  ± 1  (l>f,rmr  in  units  of  pixels.  For  example, 
if  the  SST  took  an  image  through  a  turbulent  atmosphere  with  a  r0  value  of  10cm,  using  the 
model  without  Poisson  noise,  a  potential  error  from  its  true  position  in  pixel  position  would 
be  ±  1.0917  pixels.  Assuming  the  object  being  observed  by  SST  is  in  geosynchronous  orbit 
at  approximately  36,000km,  this  would  correspond  to  about  168.44m  observed  difference 
from  its  true  location.  Note,  Eqs.  (3.56)  and  (3.57)  are  specific  for  images  observed  by  the 
SST  with  an  integration  time  of  100ms. 
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It  was  also  observed  from  the  data  that  lower  wind  speeds  (near  0 m/s)  for  a  specific 
r0  caused  higher  observed  position  errors.  The  opposite  is  also  true  that  higher  wind 
speeds  caused  lower  observed  position  errors.  However,  as  r0  increased,  the  difference 
(or  variance)  between  high  and  low  wind  speeds  decreases.  The  worst  case  situation  would 
be  a  very  turbulent  atmosphere  (low  r0  value)  and  a  low  wind  speed.  This  makes  sense 
from  a  physical  perspective  as  large  scale  turbulent  eddies  cause  tilt  [5];  if  the  wind  is  not 
blowing,  the  tilt  caused  by  that  turbulent  eddy  will  remain  in  the  field  of  view  causing  the 
observed  image  to  be  translated. 

Eqs.  (3.56)  and  (3.57)  can  be  considered  the  observation  uncertainty  and  take  into 
account  both  the  imperfections  in  the  optics  of  a  system  and  the  tilt  phase  error  due  to 
atmospheric  turbulence.  This  model  allows  for  the  observation  uncertainty,  referred  to  as 
Oq,  in  the  next  subsection,  to  change  as  a  function  of  r0. 

3.7.2  The  Total  Observation  Uncertainty. 

The  observation  uncertainty  affects  the  ability  to  accurately  determine  the  position 
of  a  detected  space  object.  The  observation  uncertainty  affects  the  ability  to  accurately 
determine  the  angular  position,  rate,  and  acceleration.  From  a  purely  deterministic 
perspective,  once  the  angular  rates  and  accelerations  are  measured  for  both  the  right 
ascension,  a,  and  declination,  6,  angles,  it  is  possible  to  determine  the  range,  p,  and  range 
rate,  p  [11].  The  two  astronomical  coordinates,  a  and  6,  combine  to  specify  the  a  point 
on  the  celestial  sphere  in  the  equatorial  coordinate  system.  Theoretically,  once  all  the 
variables  of  Eq.  (3.58)  are  measured  from  ground  based  optical  sensors  or  determined  from 
the  equations  of  motion,  then  all  variables  of  the  state  equation  of  the  orbiting  object  are 
known  and  it  is  possible  to  calculate  the  orbit  of  an  object.  The  ability  to  determine  range 
and  range  rate  is  derived  from  the  equations  of  motion  in  which  the  acceleration  of  either 
angle  is  a  mathematical  function,  where  the  function  is  represented  by  fa  and  fs,  of  the 
angular  position,  angular  rate,  range  and  range  rate.  In  Eq.  (3.58),  the  dot  and  double  dot 
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above  the  orbital  parameter  represent  the  first  and  second  derivatives  of  the  right  ascension, 
a,  declination,  6,  and  range,  p.  From  an  optical  detection  perspective,  it  is  very  difficult  to 
determine  the  range  and  range  rate  parameters;  however,  if  acceleration  can  be  measured, 
then  through  the  equations  of  motion,  the  range  and  range  rate  data  can  be  determined. 


a  =  fa(a,  a,  6,  S,p,p) 


6  =  fs(a,  a,  6, 6,p,p) 


(3.58) 


However,  from  a  practical  perspective,  there  is  a  large  uncertainty  in  measuring  the 
acceleration,  cr&0 ,  of  an  orbital  object  which  translates  into  a  large  uncertainty  in  estimating 
the  range  and  range  rate  [11].  By  using  the  equations  described  in  the  Maruskin,  et  al. 
paper,  it  is  possible  to  determine  the  error  associated  with  measuring  angular  position,  rate, 
and  acceleration  using  the  SST. 

It  is  important  to  consider  the  errors  or  observation  uncertainties  in  measuring  the 
angular  position,  rate,  and  acceleration  of  a  detected  space  object.  Assume  that  the  angular 
motion  can  be  modeled  kinematically  as  listed  in  Eq.  (3.59)  [11]. 

a{t)  =  a0  +  a0(t  -  ta )  +  -aQ{t  -  t0f  (3.59) 


Provided  there  are  N  equally  spaced  observations  (minimum  of  N  >  3  observations)  of 
the  same  object  over  a  span  of  time  T ,  it  is  possible  to  determine  the  angle  at  epoch,  a0, 
the  angular  rate  at  epoch,  a0,  and  angular  acceleration  at  epoch,  a0.  Assuming  that  the  N 
angular  measurements  are  uncorrelated  with  each  other,  a  least-squares  estimation  problem 
is  used  to  minimize  the  cost  function 

1  N 

7 = ^  Yi[a{ti)  -  ai]2>  (3-60) 

(=i 

where  ffi„  is  the  observation  uncertainty,  a,  are  the  actual  measurements,  and  ait,)  are  the 
predicted  measurements;  these  quantities  are  used  to  estimate  a0,  a0,  and  a0.  Assuming 
that  the  measurements  are  taken  at  equal  times,  then  tj  =  ta  +  (T /2 ri)i,  where  i  = 
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-77,  — (t7  —  1), -1,0, 1, 77,  creating  a  total  of  N  =  2/7  +  1  measurements  over  the  total  time 
span  T  [11].  Due  to  taking  measurements  at  equal  times,  the  odd  terms  of  the  information 
matrix  A  sum  to  zero: 


A  = 
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6  1 2n  / 

The  inverse  of  this  matrix  is  covariance  matrix  Q  which  provides  the  accuracy  of  the  angle 
quantities  measured: 


Q 


O: 


(277  +  1) 


0  15/2^ 


-¥(?) 
3(f)2  0 


2  \4 


fijr  0  45(f) 


(3.62) 


By  taking  the  square  root  of  the  diagonal  entries,  it  is  possible  to  determine  the  standard 
deviation,  if/,  of  the  estimated  angular  motion  variables  a0,  a0,  and  a0  listed  in  Eq.  (3.59): 
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(3.65) 


The  purpose  of  using  position  algorithms  that  improve  the  position  accuracy  of  a 
detected  space  object  such,  as  the  CoS  method,  is  to  reduce  the  observation  uncertainty. 
This  is  in  turn  reduces  the  errors  associated  with  making  measurements  of  angular  position, 
rate,  and  acceleration.  As  the  SST  determines  the  track  for  a  detected  object  based  on 
three  subsequent  detections,  reducing  the  observation  uncertainty  as  much  as  possible  is 
important. 
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IV.  Results 


The  results  on  this  research  are  based  on  a  statistical  analysis  of  the  calculated  angular 
rate  of  30  different  stars  using  the  SST  data.  The  data  is  captured  by  the  SST  over  three 
days;  10  random  stars  are  used  from  each  day  of  data.  As  the  purpose  of  this  research 
is  to  determine  which  position  algorithm  returns  the  most  accurate  angular  rate,  the  stars 
selected  are  bright;  thus,  the  ability  to  detect  a  space  object  has  no  impact  in  determining 
the  position.  In  this  research,  on  average,  18.9  frames  of  sequential  data  is  taken  for  each 
star  in  order  that  a  mean  with  some  statistical  significance  could  be  determined.  Lastly, 
the  window  is  limited  to  only  one  star;  this  is  to  avoid  issues  where  multiple  objects  are 
detected  in  a  single  frame  of  data,  which  would  manipulate  the  centroid  operations.  In 
full  frame  data,  an  unwindowed  frame  of  data  may  have  multiple  detected  objects.  The 
ability  to  determine  the  angular  rate  for  multiple  detected  objects  would  involve  a  clustering 
algorithm;  however,  this  research  is  limited  to  measuring  the  angular  rate  of  a  single  space 
object. 

4.1  Assessing  the  Position  and  Angular  Rate  Algorithms 

The  results  of  this  research  show  that  from  a  statistical  perspective,  that  all  angular  rate 
algorithms,  given  18.9  frames  of  data  on  average,  will  calculate  a  mean  angular  rate  ratio, 
E{ji\,  for  the  all  30  stars  close  to  the  actual  angular  rate  of  a  star.  As  is  noted  in  section  3.6, 
the  calculated  angular  rate  is  divided  by  the  actual  angular  rate  of  a  star;  thus,  a  perfectly 
calculated  angular  rate  divided  by  the  actual  angular  rate  would  produce  an  angular  rate 
ratio  of  1. 

Looking  at  table  (4.1)  and  considering  only  the  statistical  mean,  E\ju\,  of  the  average 
angular  rate,  /j ,  for  each  star,  the  centroid  method  performs  the  best  (1.0056864)  given 
about  19  frames  of  data  followed  in  order  by  the  Col  method  (1.0057782),  the  CoSi  method 
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Table  4.1:  Calculated  Mean  Angular  Rate  for  Each  Angular  Rate  Algorithm 


data  day 

#of  frames 

fj  centroid 

/r  Col 

yuCoSi 

A*  CoS2 

73 

23 

1.0049711 

1.0054758 

1.0056160 

1.0052673 

73 

21 

1.0070785 

1.0062551 

1.0061252 

1.0064317 

73 

13 

1.0052292 

1.0055392 

1.0042487 

1.0040030 

73 

22 

1.0059939 

1.0059618 

1.0057218 

1.0060486 

73 

20 

1.0083811 

1.0075002 

1.0074646 

1.0074760 

74 

20 

1.0053867 

1.0057544 

1.0060864 

1.0060435 

74 

18 

1.0068006 

1.0070948 

1.0072384 

1.0071367 

74 

19 

1.0067455 

1.0067682 

1.0065905 

1.0065977 

74 

17 

1.0049657 

1.0052306 

1.0064570 

1.0066316 

74 

18 

1.0049777 

1.0057517 

1.0057167 

1.0058569 

74 

14 

1.0067543 

1.0067997 

1.0066813 

1.0071213 

74 

19 

1.0057002 

1.0058825 

1.0059335 

1.0060599 

74 

20 

1.0066515 

1.0056373 

1.0049228 

1.0052492 

74 

21 

1.0073051 

1.0062502 

1.0054367 

1.0053855 

75 

18 

1.0043525 

1.0045270 

1.0047913 

1.0047076 

75 

19 

1.0055581 

1.0055253 

1.0053320 

1.0052745 

75 

19 

1.0043362 

1.0050875 

1.0048318 

1.0047935 

75 

18 

1.0044765 

1.0044657 

1.0048413 

1.0045600 

75 

19 

1.0055261 

1.0060162 

1.0062280 

1.0061452 

75 

19 

1.0076103 

1.0055851 

1.0052363 

1.0050916 

75 

17 

1.0026446 

1.0042524 

1.0056341 

1.0056372 

75 

13 

1.0034606 

1.0050959 

1.0059679 

1.0060691 

75 

19 

1.0056633 

1.0060919 

1.0061616 

1.0062177 

75 

19 

1.0055581 

1.0055253 

1.0053320 

1.0052745 

74 

21 

1.0051755 

1.0050625 

1.0051932 

1.0051160 

73 

19 

1.0062770 

1.0060259 

1.0060197 

1.0059555 

73 

19 

1.0042889 

1.0051152 

1.0055301 

1.0056028 

73 

21 

1.0052778 

1.0053031 

1.0053482 

1.0051598 

73 

21 

1.0070172 

1.0060705 

1.0058571 

1.0061884 

73 

21 

1.0064290 

1.0076935 

1.0078084 

1.0076327 

1.0056864 

1.0057782 

1.0058118 

1.0058245 

0.0012631 

0.0008212 

0.0008086 

0.0008684 
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without  noise  outlier  removal  (1.0058118),  and  lastly  by  the  CoS2  method  with  the  noise 
outlier  removal  (1.0058245).  The  performance  of  the  centroid  translates  to  an  angular  rate 
of  73.3358 lyurad/sec  compared  to  the  performance  of  the  CoS2,  which  corresponds  to  an 
angular  rate  of  7 3.3458 9/rrad/sec .  This  is  only  a  difference  of  0.0 100 8/7  rad/s  cc,  meaning 
that  all  the  methods  are  very  close.  The  statistical  mean,  E[p],  is  a  measurement  of  the 
accuracy  of  the  angular  rate  methods.  However,  looking  at  the  standard  deviation  of  the 
accuracy,  -  7[//])2],  the  centroid  method  has  the  largest  variability  at  0.0012631, 

while  the  CoSi  is  the  most  consistent  at  0.0008086.  With  that  said,  each  angular  rate 
algorithms  would  fit  well  within  1  standard  deviation  of  the  other  methods  standard 
deviation;  therefore,  the  difference  in  the  accuracy  of  each  measurement  could  be  explained 
by  the  statistical  randomness  of  the  samples  chosen  in  this  research.  Analyzing  the 

variance,  E[cr2]  of  the  calculated  angular  rates,  which  is  a  measurement  of  the  precision  of 
the  angular  rate  algorithms,  the  most  precise  angular  rate  algorithm  is  the  CoSi  algorithm 
with  an  average  variance  of  E[cr2]  -  0.0001185.  This  is  approximately  4  times  more 
precise  than  the  centroid  algorithm  at  E[cr2]  =  0.0004462.  Adding  one  standard  deviation 
to  the  mean  of  the  centroid  method  results  in  an  angular  rate  ratio  of  1.02681.  This 
corresponds  to  an  angular  rate  of  74.876  l//rad/scc.  Performing  this  same  operation  using 
the  CoS2  returns  an  angular  rate  ratio  of  1.01670,  which  corresponds  to  an  angular  rate  of 
74. 1 392/irad/scc.  The  difference  between  the  two  is  0.7369yurad/sec. 

One  method  to  quantify  the  quality  of  the  measurement  is  the  variance  of  estimation 
error  [16].  To  compare  the  algorithms  from  both  an  accuracy  and  precision  perspective, 
Eq.  (4.1)  is  used  to  quantify  the  total  error,  etotai,  where  E[cr 2]  is  the  variance,  E[/j]  is  the 
mean,  and  T  is  the  true  angular  rate  ratio.  In  this  research,  7=1.  The  total  error  will  have 
units  of  error  squared. 

e2total  =  E[<t2]  +  (EitA  -  T  f  (4.1) 
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Table  4.2:  Calculated  ±  3  Standard  Deviation  of  the  Angular  Rate  for  Each  Angular  Rate 
Algorithm 


data  day  #of  frames 

3cr  centroid 

3cr  Col 

3cr  CoS, 

3<x  CoS2 

73 

23 

0.0445912 

0.0298581 

0.0251999 

0.0289219 

73 

21 

0.0392077 

0.0296648 

0.0271911 

0.0292341 

73 

13 

0.0963049 

0.0446007 

0.0315409 

0.0354693 

73 

22 

0.0717623 

0.0382374 

0.0381151 

0.0280439 

73 

20 

0.0724458 

0.0364097 

0.0303710 

0.0280032 

74 

20 

0.0428336 

0.0334345 

0.0314792 

0.0336935 

74 

18 

0.0571243 

0.0451754 

0.0343813 

0.0383804 

74 

19 

0.0569313 

0.0316126 

0.0299362 

0.0202356 

74 

17 

0.1246202 

0.0865768 

0.0360987 

0.0269769 

74 

18 

0.0522510 

0.0271277 

0.0238602 

0.0230673 

74 

14 

0.0361371 

0.0239891 

0.0178208 

0.0195015 

74 

19 

0.0501664 

0.0336210 

0.0227849 

0.0178829 

74 

20 

0.0599039 

0.0340425 

0.0451242 

0.0433622 

74 

21 

0.0671119 

0.0452155 

0.0362285 

0.0329722 

75 

18 

0.0577331 

0.0441061 

0.0359139 

0.0433578 

75 

19 

0.0621439 

0.0335418 

0.0263784 

0.0306161 

75 

19 

0.0526860 

0.0375872 

0.0382641 

0.0367432 

75 

18 

0.0635416 

0.0325461 

0.0318498 

0.0320464 

75 

19 

0.0432926 

0.0265838 

0.0234228 

0.0280347 

75 

19 

0.0484179 

0.0317165 

0.0362577 

0.0360900 

75 

17 

0.1195076 

0.0669008 

0.0248923 

0.0267181 

75 

13 

0.0637285 

0.0488859 

0.0378567 

0.0362041 

75 

19 

0.0692308 

0.0441056 

0.0389342 

0.0328561 

75 

19 

0.0621439 

0.0335418 

0.0263784 

0.0306161 

74 

21 

0.0726213 

0.0455256 

0.0386212 

0.0383896 

73 

19 

0.0625248 

0.0342114 

0.0345918 

0.0339022 

73 

19 

0.0636214 

0.0404866 

0.0366572 

0.0408823 

73 

21 

0.0461279 

0.0334492 

0.0317194 

0.0337131 

73 

21 

0.0508438 

0.0377533 

0.0303708 

0.0322880 

73 

21 

0.0914496 

0.0603033 

0.0574650 

0.0608769 

£[3cr] 

0.0633669 

0.0396937 

0.0326569 

0.0326360 

£[rr2] 

0.0004462 

0.0001751 

0.0001185 

0.0001183 

e2 

fota/ 

0.0004785 

0.0002085 

0.0001523 

0.0001523 

&  total 

0.0218743 

0.0144379 
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0.0123399 

0.0123398 

This  method  for  describing  which  angular  rate  algorithm  is  the  most  accurate  and  precise 
from  an  equal  weighting.  This  is  method  for  quantifying  which  angular  rate  algorithm  is 
best  is  based  on  the  etotai,  which  assigns  equal  weighting  to  the  mean,  E[p],  differenced 
from  the  true  angular  rate,  T,  as  well  as  the  error  due  to  the  variance,  E[cr 2],  of  the 
calculated  angular  rate.  As  this  method  judges  the  accuracy  and  precision  from  an  equal 
footing,  the  total  error  will  tell  which  error  source  between  accuracy  and  precision  has  more 
significance  in  the  overall  calculated  angular  rate.  Note  that  the  data  in  this  research  are 
based  on  the  ratio  of  calculated  angular  rate  compared  to  the  actual  angular  rate;  a  perfectly 
calculated  angular  rate  would  yield  to  an  angular  rate  ratio  of  1.  Additionally,  the  total  error 
is  computed  in  percentage  squared.  Therefore,  the  total  error  in  these  results  is  a  percentage 
squared  from  1 . 

Based  on  this,  the  C0S1  and  CoS2  algorithms  both  have  a  total  error  of  0.01234.  This 
is  slightly  better  than  the  Col  algorithm  which  has  a  total  error  of  0.01444.  The  centroid 
method  scores  a  0.02187,  which  is  approximately  1.77  times  worse  than  the  CoS  methods. 
Given  equal  weighting  of  the  variance  and  mean,  the  precision,  E[cr 2],  is  a  more  significant 
factor  in  calculating  the  angular  rate  of  a  detected  space  object. 

Looking  at  the  total  error  for  the  centroid  algorithm  in  terms  of  angular  rate,  an  etotai 
of  0.02187  would  correspond  to  an  error  of  1 .595//rad/sec.  Looking  at  the  total  error  for 
the  CoS  methods  in  terms  of  angular  rate,  which  both  have  an  etotai  of  0.01234,  they  would 
correspond  to  an  0.899//rad/sec  in  error. 
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V.  Discussion 


The  results  of  this  research  show  that  both  the  Center  of  SNR  (CoS)  algorithms,  which 
are  based  on  the  multi-hypothesis  test  (MHT)  detection  algorithms  with  and  without  noise 
outlier  removal  performed  the  best  in  terms  of  precision  and  total  error  when  measuring  the 
angular  rate  of  an  object.  The  centroid  method  is  the  most  accurate;  however,  its  slightly 
better  performance  in  accuracy  is  well  within  one  standard  deviation  of  the  data.  This 
means  that  this  slightly  better  performance  in  measuring  the  mean  angular  rate  could  be 
explained  by  variances  in  the  data;  therefore,  a  definitive  conclusion  about  it  being  more 
accurate  angular  rate  algorithm  cannot  be  made  in  this  research. 

5.1  Relevance  of  Results 

The  relevance  of  these  results  are  that  the  binary  hypothesis  and  centroid  method 
currently  employed  on  the  SST  for  detecting  space  objects  and  determining  their  position 
and  angular  rate  information  could  be  improved  upon  by  implementing  the  MHT  detection 
algorithm  and  one  of  the  CoS  algorithm  for  determining  position  and  angular  rate  described 
in  this  research.  Looking  at  the  total  error  as  described  in  section  (4.1),  both  CoS  algorithms 
had  a  total  error  of  0.01234,  while  the  method  currently  employed  on  the  SST  has  a  total 
error  of  0.02187.  This  means  that  in  terms  of  total  error,  that  the  CoS  algorithms  are 
approximately  1.77  times  better  than  the  currently  employed  centroid  method.  The  total 
error  in  terms  of  angular  rate  for  the  centroid  method  is  1 .595//rad/scc.  In  comparison, 
the  CoS  algorithms  have  a  total  error  of  0.899//rad/sec.  As  the  SST  uses  around  3  to  5 
sequential  observations  to  determine  a  track  for  a  detected  object,  having  the  most  precise 
and  accurate  angular  rate  measurement  is  critical. 
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5.2  Future  Work 


The  purpose  of  this  research  is  in  categorizing  which  angular  rate  algorithm  is  accurate 
and  precise.  The  next  step  in  this  research  would  be  to  use  the  angular  rate  measurements 
from  a  detected  space  object  to  determine  its  track  and  orbital  parameters.  Additionally,  if 
this  research  is  to  be  implemented  on  the  SST,  the  code  will  need  to  be  streamlined  in  order 
that  it  runs  near  real-time.  This  requires  understanding  the  full  SST  system  in  order  that 
the  MHT  detection  and  CoS  angular  rate  algorithms  can  be  implemented  into  the  current 
operational  system.  Finally,  more  work  is  required  in  being  able  to  implement  this  code 
with  a  full  SST  data  frame  which  can  contain  multiple  objects  each  moving  at  a  unique 
angular  rate. 

5.3  Conclusion 

Both  CoS  methods  performed  in  this  research  are  better  at  determining  the  position 
and  angular  rate  of  a  detected  object  than  the  currently  employed  centroid  method  on  the 
SST.  This  research  demonstrated  that  improvements  can  be  made  by  using  the  MHT 
methods.  The  implications  are  that  the  observation  uncertainty  will  go  down,  thereby 
improving  the  estimated  angular  position  and  rate.  This  in-turn  improves  the  ability  to 
determine  the  track  of  a  detected  space  object  and  improve  the  ability  to  perform  SSA. 
Additional  research  will  need  to  be  conducted  to  quantify  by  how  much  tracks  of  detected 
space  objects  are  improved  by  when  using  the  MHT  algorithms  in  conjunction  with  the 
CoS  algorithms. 
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Appendix:  MATLAB  code 


A.l  Main  MATLAB  code 
%Anthony  SI i gar 
%SST  detection 


clear ; 

close  all; 

%Day 

seeing 

%73 

5.9 

%74 

6.1 

%75 

5.6 

%81 

3.5 

%82 

4.8 

%83 

6.4 

clc ; 


(cm) 


Bx=2 . ®2843618486®734e+®3 ;  %x_correct=x_window+view_xl+B 

By=3 .0753345 1883 3674e+®3 ;  %y 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


load  data_cube073_3d;  %data  used 

rdata=data_cube073_3d; 

load  hypo59;  %loads  hypotheses  with  5.9  cm  seeing  parameter 


%starl 

framel=552;  %first  frame  to  view 
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frame2=575 ; 

%last  frame 

to 

view_xl=60 ; 

%coordinate 

of 

view_x2=80 ; 

%coordinate 

of 

view_yl=l ; 

%coordinate 

of 

view_y2=20© ; 

%coordinate 

of 

view 

1st  x  value  in  viewing  range 
last  x  value  in  viewing  range 
1st  y  value  in  viewing  range 
last  y  value  in  viewing  range 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%sets  size  of  psf  to  the  size  of  the  window 
if  (view_x2-view_xl)>=20 

psf=padarray(hypo , [view_y2-view_yl-19  view_x2-view_xl-19] ,  0, ’pre’); 
z=psf ( : , : , 5) ; 

[num  indx]=max(z(:)) ; 

[Yshift  Xshift]  =  ind2sub(size(z) , indx) ;  %uses  index  value  to  find  location  of 
psf=circshift(psf , [1-Yshift, 1-Xshift]) ; 
else  psf=padarray(hypo , [view_y2-view_yl- 19  ©],©,’pre’); 
psf=psf ( : , 1 : (view_x2-view_xl+l) , : ) ; 
z=psf ( : , : , 5) ; 

[num  indx] =max(z( : )) ; 

[Yshift  Xshift]  =  ind2sub(size(z) , indx) ;  %uses  index  value  to  find  location  of 
psf=circshift(psf , [1-Yshift, 1-Xshift]) ; 

end 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 
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%detection  tests  and  position  array 
for  ii=framel : frame2 

Dsst=rdata(view_yl : view_y2 , view_xl : view_x2 , ii) ; 
[xcent(ii)  ycent(ii)  HQ]=centroid(Dsst) ; 
[xCOI(ii)  yCOI(ii)  Hl]=COI(Dsst) ; 

[xMHT(ii)  yMHT(ii)  H3]=MHT3(Dsst,  psf) ; 
[xSNR(ii)  ySNR(ii)  H4]=MHT4(Dsst,  psf); 

end 


%load  SST  frame 
%centroid  function  call 
%CoI  function  call 
%CoSl  function  call 
%CoS2  function  call 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%position  information 
Xlocl=xcent+view_xl+Bx ; 
Ylocl=ycent+view_yl+By ; 


%centroid  method 
%centroid  method 


XlocCOI=xCOI+view_xl+Bx ; 
YlocCOI=yCOI+view_y 1+By ; 


%CoI  method 
%CoI  method 


XlocSNR=xSNR+view_xl+Bx ; 
YlocSNR=ySNR+view_y 1+By ; 


%CoS2  method 
%CoS2  method 


XlocMHT=xMHT+view_xl+Bx ; 
YlocMHT=yMHT+view_y 1+By ; 


%CoSl  method 
%CoSl  method 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%converts  pixels  to  radians  in  field  of  view 
radperpix=tan((30/3 .49)*10''-6)  ;  %SST  rad/pixel 
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sid_radsec=(2*pi/(23 . 9344699*3600)) ;  %sidereal  rate  (rad/sec) 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


for  kk=framel : frame2-l 

xslope (kk) = (Xloc 1 (kk+1) -Xloc 1 (kk) ) *radperpix ; 
yslope (kk) = (Yloc 1 (kk+1) -Yloc 1 (kk) ) *radperpix ; 


%centroid  method 
%centroid  method 


xslopeCOI(kk)=(XlocCOI(kk+l) -XlocCOI(kk)) *radperpix ;  %CoI  method 
yslopeCOI(kk)=(YlocCOI(kk+l) -YlocCOI(kk)) *radperpix ;  %CoI  method 


xslopeSNR(kk)=(XlocSNR(kk+l) -XlocSNR(kk)) *radperpix ;  %CoS2  method 
yslopeSNR(kk)=(YlocSNR(kk+l) -YlocSNR(kk)) *radperpix ;  %CoS2  method 


xslopeMHT (kk)=(XlocMHT (kk+1) -XlocMHT (kk) ) *radperpix ;  %CoSl  method 
yslopeMHT (kk)=(YlocMHT (kk+1) -YlocMHT (kk) ) *radperpix ;  %CoSl  method 


end 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%centroid  method 

xslopel=xslope . /sid_radsec ;  %ratio  between  object’s  angular  rate  to  sidereal  rate 
yslopel=yslope . /sid_radsec ; 

combinedmag=(sqrt ((xslope . ~2)+ (yslope . "2))) . /sid_radsec ; 
avgCen=ones(l , length(combinedmag))*mean(combinedmag(combinedmag~=Q)) ; 
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pstdCen=avgCen+(ones(l ,  length(combinedmag))*3*std(combinedmag(combinedmag~=0)))  ; 
nstdCen=avgCen-  (ones(l ,  length(combinedmag))*3*std(combinedmag(coinbinedmag~=0)))  ; 

%CoI  method 

xslopeCOIl=xslopeCOI . /sid_radsec ;  %ratio  between  object’s  angular  rate  to  sidereal 
yslopeCOI l=yslopeCOI . /sid_radsec ; 

combinedmagCOI=(sqrt((xslopeCOI . ~ 2) + (yslopeCOI . "2))) . /sid_radsec ; 
avgCOI=ones ( 1 ,  length (combinecLmag)  )  *mean(combinedmagCOI  (combinecLmagCOI  ~=0)  )  ; 
pstdCOI=avgCOI+(ones(l , length(combinedmag))  *3*std(combinedmagCOI(combinedmagCOI~=0) 
nstdCOI=avgCOI- (ones(l , length(combinedmag))  *3*std(combinedmagCOI(combinedmagCOI~=0) 

%CoS2  method 

xslopeSNRl=xslopeSNR./sid_radsec;  %ratio  between  object’s  angular  rate  to  sidereal 
yslopeSNRl=yslopeSNR . /sid_radsec ; 

combinedmagSNR=(sqrt((xslopeSNR.  ,'2)  +  CyslopeSNR.  ''2)))  . /sid_radsec ; 
avgSNR=ones(l , length(combinedmag))*mean(combinedmagSNR(combinedmagSNR~=0)) ; 
pstdSNR=avgSNR+(ones(l , length(combinedmag))*3*std(combinedmagSNR(combinedmagSNR~=0) 
nstdSNR=avgSNR- (ones(l , length(combinedmag))*3*std(combinedmagSNR(combinedmagSNR~=0) 

%CoSl  method 

xslopeMHTll=xslopeMHT./sid_radsec;  %ratio  between  object’s  angular  rate  to  sidereal 
yslopeMHTl l=yslopeMHT . /sid_radsec ; 

combinedmagMHT=(sqrt((xslopeMHT.  ,'2)  +  CyslopeMHT.  ''2)))  . /sid_radsec ; 
avgMHT=ones ( 1 ,  length(combinedmag)  )  ''meanCcombinedmagMHT (combinedmagMHT~=0)  )  ; 
pstdMHT=avgMHT+(ones(l , length(combinedmag))  "'3*std(combinedmagMHT (combinedmagMHT~=0) 
nstdMHT=avgMHT- (ones(l , length(combinedmag))  "'3*std(combinedmagMHT (combinedmagMHT~=0) 
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o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%plotting 

figure(l) 

subplot(3 ,2,1) 

set(gca,  ’FontSize’,  12); 

plot (Xlocl ,’ color ’ ,  ’b ’ , ’ line style ;hold  on; plot(XlocCOI , ’ color ’ , ’m’ , ’ linesty 

xlim([  framel  frame2  ]) 

xlabel( ’ frame  number’); 

ylabelC ’pixel  x-direction’ ) ; 

title( ’ Position  (x)’) 

legend ( ’ Centroid’ , ’ Col ’ , ’CoS’ , ’CoS  w/outlier  removal ’ , ’ location’ , ’ southeast ’ ) ; 

subplot(3 ,2,2) 
set(gca,  ’FontSize’,  12); 

plot (Ylocl ,’ color ’ ,  ’b ’ , ’ line style ’,’--’) ;hold  on; plot(YlocCOI , ’ color ’ , ’m’ , ’ linesty 

xlim([  framel  frame2  ]) 

xlabel( ’ frame  number’); 

ylabelC ’pixel  y-direction’ ) ; 

title( ’ Position  (y)’) 

legendC ’ Centroid’ ,’ Center  of  Intensity ’, ’MHT’ , ’MHT  w/outlier  removal ’,’ location’ ,’ S' 

subplot(3 ,2,3) 
set(gca,  ’FontSize’,  12); 

plot (xslope , ’ color’ ,  ’ b ’ , ’ line style ’ , ’ -- ’ ) ;hold  on; plot (xslopeCOI , ’ color ’ , ’m’ , ’ line 
xlim([  framel  frame2  ]) 
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xlabelC ’ frame  number’); 
ylabelC’rate  (radians/frame)  (x)’); 
title( ’Angular  rate  radians/frame  (x)’) 

legend ( ’ Centroid’ , ’ Col ’ , ’CoS’ , ’CoS  w/outlier  removal ’ , ’ location’ , ’ southeast ’ ) ; 

subplot(3 ,2,4) 
setCgca,  ’FontSize’,  12); 

plot(yslope , ’ color’ ,  ’ b ’ , ’ line style ’ , ’ -- ’ ) ;hold  on; plot (yslopeCOI , ’ color ’ , ’m’ , ’ line 

xlim([  framel  frame2  ]) 

xlabelC ’ frame  number’); 

ylabelC’rate  Cradians/frame)  Cy)’); 

titleC ’Angular  rate  radians/frame  Cy)’) 

legend C ’ Centroid’ , ’ Col ’ , ’CoS’ , ’CoS  w/outlier  removal ’ , ’ location’ , ’ southeast ’ ) ; 

subplotC3 ,2,5) 
setCgca,  ’FontSize’,  12); 

plotCxslopel , ’ color ’ ,  ’ b ’ , ’ line style ’,’--’) ;hold  on; plot CxslopeCOIl , ’ color ’ , ’m’ , ’ li 
xlimCL  framel  frame2  ]) 
xlabelC ’ frame  number’); 

ylabelC’ ratio  rad/frame/sidereal  rate  C x)’); 
titleC ’ radians/frame  divided  by  sidereal  rate  Cx)’) 

legend C ’ Centroid’ , ’ Col ’ , ’CoS’ , ’CoS  w/outlier  removal ’ , ’ location’ , ’ southeast ’ ) ; 

subplotC3 ,2,6) 
setCgca,  ’FontSize’,  12); 

plot Cyslopel , ’ color ’ ,  ’ b ’ , ’ line style ’,’--’) ;hold  on; plot CyslopeCOIl , ’ color ’ , ’m’ , ’ li 
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xlim([  framel  frame2  ]) 
xlabelC ’ frame  number’); 

ylabelC’ ratio  rad/frame/sidereal  rate  (y)’); 
title( ’ radians/frame  divided  by  sidereal  rate  (y )’) 

legend ( ’ Centroid’ , ’ Col ’ , ’CoS’ , ’CoS  w/outlier  removal ’ , ’ location’ , ’ southeast ’ ) ; 

figure(2) 

set(gca,  ’FontSize’,  12); 

plot (combinedmag , ’ color ’ ,  ’b ’ , ’ line style ’ , ’ -- ’ ) ;hold  on; plot(combinedmagCOI , ’ color ’ 
xlim([  framel  frame2  ]) 
xlabelC ’ frame  number’); 

ylabelC ’Magnitude  of  object  in  Cradians/frame)/Csidereal  rate)’); 
titleC ’Magnitude  of  angular  rate  divided  by  sidereal  rate’) 

legend C ’ Centroid’ , ’ Col ’ , ’CoS’ , ’CoS  w/outlier  removal ’ , ’ location’ , ’ southeast ’ ) ; 

figureC3) 

setCgca,  ’FontSize’,  12); 

plot CavgCen, ’b ’ , ’ line style ’ , ’ -- ’ ) ;hold  on;plotCavgCOI , ’m’ , ’ line style ’,’:’) ;hold  on; 
plot CpstdCen, ’b ’ , ’ line style ’,’--’) ;hold  on; plotCnstdCen, ’b ’ , ’ line style ’,’--’) ;hold  i 
plotCpstdCOI , ’m’ , ’ linestyle ’,’:’) ;hold  on;plotCnstdCOI , ’m’ , ’ linestyle ’,’:’) ;hold  on 
plotCpstdMHT, ’r’ , ’linestyle’ , ’) ;hold  on; plotCnstdMHT , ’r’ , ’linestyle’ , ’) ;hold  i 
plotCpstdSNR, ’k’) ;hold  on;plotCnstdSNR, ’k’) ;hold  on; 
xlimCL  framel  frame2  ]) 
xlabelC ’ frame  number’); 

ylabelC ’Magnitude  of  object  in  radians/frame’); 

titleC ’Magnitude  of  angular  rate  divided  by  sidereal  rate  using  mean  and  stdev’); 
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legend ( ’ Centroid’ , ’ Col ’ , ’CoS’ , ’CoS  w/outlier  removal ’ , ’ location’ , ’ southeast ’ ) ; 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


format  long ; 

ZZ=[max(avgCen)  max(avgCOI)  max(avgMHT)  max(avgSNR)  max(pstdCen)  max(pstdCOI)  max(p 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


A.2  MATLAB  code  functions:  Centroid 
function  [xpos  ypos  HI]  =  centroid(Dsst) ; 


tau=6;  %threshold  value 


0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/ 0/0/0/ 
/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%point  detection  test 

Bh=median(Dsst( : )) ;  %this  is  the  noise  determined  from  Dsst 
sigma=sqrt(mean(mean((Dsst-Bh) . "2))) ;  %takes  the  value  of  all  the  pixels 
LRT=((Dsst-Bh)/sigma) ;  %LRT  is  likelihood  ratio  test 


[m,n]=size(Dsst) ;  %m  is  size  in  y  direction,  n  is  size  in  x  direction 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%create  matrix  of  detection  points 
for  m=l:m 
for  n=l:n 


if  LRT(m,n)>tau 

Hl(m,n)=l;  %sets  detection  parameter  HI  to  1  if  there  is  an  object 
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else  Hl(m,n)=Q; 
end 


end 

end 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%centroid  operation 
xsum=sum(Hl) ; 
ysum=sum(Hl ’ )  ; 
Hlsum=sum(sum(Hl)) ; 


for  k=l:n  %centroid  operation;  finds  x  position 
xprod(k)=xsum(k) *k ; 

end 


for  1=1 :m  %centroid  operation;  finds  y  position 
yprod(l)=ysum(l) *1 ; 

end 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%returns  position  coordinates 
xpos=sum(xprod)/Hlsum;  %x 
ypos=sum(yprod)/Hlsum;  %y 


of  values 

position  of  object 
position  of  object 


end 
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A.3  MATLAB  code  functions:  Col 


function  [xpos  ypos  Hli]  =  COI(Dsst); 


tau=6;  %threshold  value 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%point  detection  test 

Bh=median(Dsst( : )) ;  %this  is  the  noise  determined  from  Dsst 
sigma=sqrt(mean(mean((Dsst-Bh) . "2))) ;  %takes  the  value  of  all  the  pixels 
LRT=((Dsst-Bh)/sigma) ;  %LRT  is  likelihood  ratio  test 


[m,n]=size(Dsst) ;  %m  is  size  in  y  direction,  n  is  size  in  x  direction 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%create  matrix  of  detection  points 
for  m=l:m 
for  n=l:n 


if  LRT(m,n)>tau 

Hl(m,n)=l;  %sets  detection  parameter  HI  to  1  if  there  is  an  object 
else  Hl(m,n)=Q; 
end 


end 

end 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%centroid  operation 
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Hli=Dsst . *H1 ;  %sets  values  in  HI  equal  to  1  back  to  intensity  values  in  Dsst 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%Center  of  Intensity  -  COI  method  1 


xsum=sum(Hli) ; 
ysum=sum(Hli ’ ) ; 
Hlsum=sum(sum(Hli)) ; 


for  k=l:n 

xprod(k)=xsum(k) *k ; 

end 

for  1=1  :m 

yprod(l)=ysum(l) *1 ; 

end 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%returns  position  coordinates 
xpos=sum(xprod)/Hlsum;  %x 
ypos=sum(yprod)/Hlsum;  %y 


of  values 

position  of  object 
position  of  object 


end 

A.4  MATLAB  code  functions:  MHT3 

function  [xpos  ypos  MHT_SNR1]  =  MHT3(Dsst,  psf) ; 
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tau=6.2212;  %threshold  value 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%LRT  detection  test 

Bh=median(Dsst( : )) ;  %this  is  the  noise  determined  from  Dsst 
sigma=sqrt(mean(mean((Dsst-Bh) . "2))) ;  %takes  the  value  of  all  the  pixels 
LRT=((Dsst-Bh)/sigma) ;  %LRT  is  likelihood  ratio  test 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%MHT 


for  kk=l:9; 


SNR(: , : , kk)=(ifft2 (fft2 (LRT) .*fft2(psf(: , : ,kk))))/sqrt(sum(sum(psf( : , : ,kk) . ~ 


end 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 

%filter  (set  to  0)  anything  that  is  less  than  tau 
[m,n, o]=size(SNR) ; 
for  c=l:o; 
for  a=l:m; 
for  b=l:n; 


if  SNR(a,b , c)>tau 

SNR(a,b,c)=SNR(a,b,c) ;  %sets  detection  parameter  HI  to  1  if  there  is  an  obj 
else  SNR(a,b,c)=0; 
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end 


end 

end 

end 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%create  MHT_SNR  combined  matrix 
MHT_SNR=zeros(m*2 ,n*2) ; 

MHT_SNR(2 : 2 :m*2 ,2:2 :n*2)=SNR( : , : , 5) ; 
MHT_SNR(2 : 2 :m*2 ,3:2 :n*2+l)=SNR( : , : , 6) ; 
MHT_SNR(3 : 2 :m*2+l ,2:2 :n*2)=SNR( : , : ,8) ; 
MHT_SNR(3:2:m*2+l,3:2:n*2+l)=SNR(: , : ,9); 
MHT_SNR=padarray (MHT_SNR , [10  10] ,0, ’both’); 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%filter  out  all  but  the  brightest  star 
[m,n]=size(MHT_SNR) ; 
MHT_SNRl=zeros(m,n) ; 
if  max (max (MHT_SNR) ) ==0 ; 
xpos=-l ; 
ypos=-l ; 

else 


[num  indx] =max (MHT_SNR ( : ) ) ; 

[Ymax  Xmax]  =  ind2sub(size(MHT_SNR) , indx) ; 
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MHT_SNRl(Ymax-5 : Ymax+5 ,Xmax-5 :Xmax+5)=MHT_SNR(Ymax-5 : Ymax+5 ,Xmax-5 :Xmax+5) ; 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%centroid  operation 
xsum=sum(MHT_SNRl)  ; 
y  sum=  sum (MHT_SNR 1 ’ ) ; 
H2sum=sum(sum(MHT_SNRl)) ; 


[y,x]=size(MHT_SNRl) ; 

for  k=l:x 

xprod(k)=xsum(k) *k ; 

end 

for  l=l:y 

yprod(l)=ysum(l) *1 ; 

end 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%return  position  values 
xpos=sum(xprod)/(2*H2sum) -5 . 5623 ; 
ypos=sum(yprod)/(2*H2sum) -5 . 7067 ; 


end 

A.5  MATLAB  code  functions:  MHT4 

function  [xpos  ypos  MHT_SNR1]  =  MHT4(Dsst,  psf) ; 
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tau=6.2212;  %threshold  value 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%LRT  detection  test 

Bh=median(Dsst( : )) ;  %this  is  the  noise  determined  from  Dsst 
sigma=sqrt(mean(mean((Dsst-Bh) . "2))) ;  %takes  the  value  of  all  the  pixels 
binmap=((Dsst-Bh)<tau*sigma) ;  %creates  a  binary  map  where  1  is  equal  to  anythin 
sigma2=sqrt  (sum(sum((binmap .  *(Dsst-Bh))  .  ''2))/sum(sum(binmap)))  ;  %creates  a  new 
LRT=((Dsst-Bh)/sigma2) ;  %LRT  is  likelihood  ratio  test 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%MHT 


for  kk=l:9; 


SNR(: , : , kk)=(ifft2 (fft2 (LRT) .*fft2(psf(: , : ,kk))))/sqrt(sum(sum(psf( : , : ,kk) . ' 


end 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 

%filter  (set  to  0)  anything  that  is  less  than  tau 
[m,n, o]=size(SNR) ; 
for  c=l:o; 
for  a=l:m; 
for  b=l:n; 
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if  SNR(a,b , c)>tau 

SNR(a,b,c)=SNR(a,b,c) ;  %sets  detection  parameter  HI  to  1  if  there  is  an  obj 

else  SNR(a,b,c)=0; 

end 

end 

end 

end 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%create  MHT_SNR  combined  matrix 
MHT_SNR=zeros(m*2 ,n*2) ; 

MHT_SNR(2 : 2 :m*2 ,2:2 :n*2)=SNR( : , : , 5) ; 
MHT_SNR(2 : 2 :m*2 ,3:2 :n*2+l)=SNR( : , : , 6) ; 
MHT_SNR(3 : 2 :m*2+l ,2:2 :n*2)=SNR( : , : ,8) ; 
MHT_SNR(3 : 2 :m*2+l ,3:2 :n*2+l)=SNR( : , : , 9) ; 
MHT_SNR=padarray (MHT_SNR , [10  10] ,0, ’both’); 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%filter  out  all  but  the  brightest  star 
[m,n]=size(MHT_SNR) ; 
MHT_SNRl=zeros(m,n) ; 
if  max (max (MHT_SNR) ) ==0 ; 
xpos=-l ; 
ypos=-l ; 

else 
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[num  indx] =max (MHT_SNR ( : ) ) ; 

[Ymax  Xmax]  =  ind2sub(size(MHT_SNR) , indx) ; 

MHT_SNRl(Ymax-5 : Ymax+5 ,Xmax-5 :Xmax+5)=MHT_SNR(Ymax-5 : Ymax+5 ,Xmax-5 :Xmax+5) ; 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%centroid  operation 
xsum=sum(MHT_SNRl)  ; 
ysum=sum(MHT_SNRl ’ ) ; 
H2sum=sum(sum(MHT_SNRl)) ; 


[y,x]=size(MHT_SNRl) ; 

for  k=l:x 

xprod(k)=xsum(k) *k ; 

end 

for  1=1 :y 

yprod(l)=ysum(l) *1 ; 

end 


o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/ 

/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o 


%return  position  values 
xpos=sum(xprod)/(2*H2sum)-5 . 5810; 
ypos=sum(yprod)/(2*H2sum)-5 . 655® ; 

end 

end 
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